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In the course of this contract a number of individual tasks have 


been performed in direct support of NASA-HSFC 1 s airborne optical com- 
munications experiment (project AVLOC) and a proposed satellite optical 
communications test. All work performed on this contract has been pre- 
viously reported either in interim technical reports or by informal re- 
ports delivered to the MSFC-COR. In this final report we summarize the 
more significant accomplishments of this project as follows: 


Appendix 1. 

Appendix 2. 

Appendix 3. 

Appendix 4. 
Appendix 5. 


Threshold Detection in an On-Off Binary Optical 
Communications System with Atmospheric Scintillation 
Near Field Antenna Patterns of Obstructed Casse- 
grainian Telescopes 

Study of Atmospheric Effects on Laser Communications 
Systems 

Requirements Study for the AVLOC Experiment 
Data Analysis for AVLOC Project 


In addition to the preceding tasks, the project director has (1) pre- 
pared a study on the state of the art in atmospheric communications, (2) 
prepared a measurements program document for the AVLOC project, (3) served 
as co-principal investigator on both the AVLOC and ATS-G communications 
experiments, and performed various other tasks at the request of the pro- 
ject contracting officer’s technical representative. 



APPENDIX 1 


THRESHOLD DETECTION IN AN 
ON-OFF BINARY COMMUNICATIONS CHANNEL 
WITH ATMOSPHERIC SCINTILLATION 


ABSTRACT 


The optimum detection threshold in an on-off binary optical commu- 
nications system operating in the presence of atmospheric turbulence has 
been investigated assuming a poisson detection process and log normal 
scintillation. The dependence of the probability of bit error on log 
amplitude variance and received signal strength has been analyzed and 
semi-emperical relationships to predict the optimum detection threshold 
derived. On the basis of this analysis a piecewise linear model for an 
adaptive threshold detection system is presented. The bit error proba- 
bilities for non-optimum threshold detection systems have also been 
investigated. 
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INTRODUCTION 


It is well known that atmospheric scintillation not only in- 
creases the bit error probability in a pulse code modulated optical 
communications channel but also influences the decision level for 
optimum threshold detection. The later effect is of considerable 
practical importance to the design of efficient pulse code modulated 
optical cormuni cations systems regardless of whether fixed or adap- 
tive threshold detection is used. Fried and Schmelzer^ have consid- 
ered the bit error rates in an optical communications channel by as- 
suming gauss i an detection statistics, an approximation which is valid 
only for large numbers of signal photons. Titterton and Speck^ have 
treated the problem using Poisson statistics so that their results 
are valid for small numbers of signal photons also. Although both 
previous investigators have recognized that the optimum detection 
threshold changes in the presence of scintillation, neither has 
studied the effect in depth. 

In this paper we investigate the optimum detection threshold in 
an on-off binary optical communications channel as a function of the 
number of received photoelectrons and the strength of scintillation. 
Poisson detection statistics and log normal scintillation are assumed 
and, on the basis of this model, expressions to predict the optimum 
detection threshold are derived. In addition, we have investigated 
the bit error probabilities for sub-optimum choices of detection 
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threshold. This is of importance to the design and analysis of 
optical communications systems since exact optimization of the 
threshold level is never possible, especially in the presence of 
atmospheric effects. Application of these results to the design of 
fixed and adaptive threshold optical receivers is discussed. 

ERROR PROBABILITY IN THE OPTICAL CHANNEL. 

We shall consider an on-off binary optical communications chan- 
nel in which the laser transmitter may send either a 'T, which cor- 
responds to a pulse being transmitted, or a 'O' which corresponds to 
no pulse transmitted. Let S be the number of signal photoelectrons 
per pulse received at the detector when a '1' is transmitted and eS 
be the number of received signal photoelectrons when a 'O' is trans- 
mitted. Here e is the reciprocal of the modulator extinction ratio. 
We will assume that in addition to the signal photoelectrons the 
detector receives N noise photoelectrons per pulse, mostly due to 
background. In a threshold detection system the receiver interprets 
the received signal as either a 'V or a '0' depending on whether the 
total number of received photoelectrons is greater or less than some 
threshold level T. 

In the presence of scintillation the probability of a detection 
error is 


P E = P{0)P FA + p O)0-P D ) 


( 1 ) 
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Where P(0) and P(l) are the apriori probabilities of sending 'O' and a 
'V respectively, Pp^ is the false alarm probability, ie. the prob- 
ability of a T being received given that a 'O' was sent, and Pp is 
the detection probability, ie. the probability of a '1' being re- 
ceived given that a /l T was sent, Pp and Pp^ are given by 


FA 


oo 

-( t. 

I j-t 


exp [-(N+eS)] 

J • 


f(S)dS 


( 2 ) 


03 

2 exp [- ( N+S ) ] f(S)dS 

Jn 3=t J - 


(3) 


Assuming log normal scintillation and a symmetrical pulse code we 
obtain 


P = - + 1 e 

E 8C„ 


C £ 2T(T-1) 






v^i */ T « wr«i<T-i/2)) 

6 T (T /dZ-j 


- T / 


- -8C„ Z 


2 

£ c 0 * 


+ ( (Sq-N)e+N) J e ‘ u Y (T,((S 0 -N)e+N)e 


(Z 0 +4C £ (T-l/2)) 


)dZ f 


( 4 ) 


where y is the incomplete gamma function 


Y* (T,n) = \ S e~ n n X ,, 

n T X=T ~JT (! 

Sq is the average number of signal photoelectrons per pulse and C £ 

is the variance of ln(S). 
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The bit error probability has been evaluated by numerical in- 
tegration of equation 4 on a UNIVAC 1108 computer. The threshold T 
was initially taken to be the no scintillation optimum threshold 
given by 3 


T„ = 


_ S ( 1-e) + In 


ln( (S+N)/{N+eS 


(P(0)(P(D) 

N+eS)) 


( 6 ) 


An iterative procedure was then used to minimize the error probability 
as a function of T, while holding Sg, N, e and constant. In this 
way we were able to simultaneously determine the optimum threshold in 
the presence of scintillation and study the effect of sub-optimum 
threshold on the error rates. Calculations were repeated for values 
of the parameters Sg, e, N and C £ over the ranges likely to be 
encountered in practice. 

Figure 1 shows the probability of a bit error at optimum thresh- 
old as a function of the number of received signal photoelectrons and 
the log amplitude variance of scintillation with one background photo- 
electron and an extinction ratio of 15 db. The computed values of 

bit error probability agree with those previously published by 

2 

Titterton and Speck . 

One feature of figure 1 which deserves note is the difference 
in the slope of the curves for large and small values of the log 
amplitude variance. When is less than about 0.02 the error prob- 
ability decreases rapidly with increasing number of signal photo- 
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electrons, whereas for large values of C the error probability tends 
to become relatively independent of the number of signal photoelectrons 
provided that Sq is not too small. Thus there appears to be two re- 
gions, one in which the error probability is determined primarily 
by scintillation and the other in which it is mainly determined by 
signal strength. This dicotomy, which we might call the signal 
limited and scintillation limited cases, will be seen more clearly 
when we consider the optimum detection threshold in the presence of 
scintillation. 

An alternate method of displaying these results is in terms of 
the transmitted power required to achieve a given bit error rate in 
the presence of scintillation as compared to the power required when 
there is no scintillation. In figure 2 the required power margin is 
shown as a function of bit error rate and log amplitude variance. 

These results are essentially independent of the choice of N and e 
as long as both are relatively small. Fried and Schmeltzer^ have 
developed a similar set of curves based on the assumption of gaussian 
detection statistics. Comparison of our results with those of Fried 
and Schmeltzer shows that the two models agree reasonably well in 
predicting the required power margin for values of C & less than about 
0.02, but that the gaussian model may underestimate the necessary 
power margin by as much as 8 db for C equal to 0.03 and 30 db for 

/V 

c 

equal to .05 at error rates of 10” . For a log amplitude variance 
of 0.1 the difference in the power margins predicted by the gaussian 
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and Poisson models may exceed 80 db. Thus the gaussian model is 
reasonably accurate in the case of weak scintillation but gives an 
extremely over-optimistic estimate of a communications system 
performance in the presence of strong scintillation. 

Since the Poisson distribution approaches a gaussian distri- 
bution as the mean increases, it might be expected that our results 
should reduce to those of Fried and Schmeltzer in the case of a large 
average number of received signal photoelectrons. However, for the 
system under consideration a large value of Sq is not sufficient to 
insure that the Poisson counting statistics can be adequately ap- 
proximated by a gaussian distribution. The means of the Poisson 
distributions in question are not S Q but (N + S) and (N + e$), where 
the variable S is averaged over all values from zero to infinity. 

Thus if N is small, the mean of the Poisson distribution will be 
small over part of the range of S, even though Sq is large. Clearly, 
when is large the contribution to the error probability from that 
part of the range of S in which the gaussian approximation is invalid 
will be greater. Physically this means that during a deep fade there 
are only a few photoelectrons per pulse; therefore gaussian statistics 
can not be used. A second, and perhaps more significant, problem 
with Fried and Schmeltzer's model is the assumption of additive noise 
independent of the signal strength. This assumption is equivalent to 
approximating the Poisson distribution by a gaussian distribution with 
constant variance. Since the variance of a Poisson distribution is 
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equal to the mean, this assumption is valid only if 

N+S-N+eS- constant (7) 

for all values of S with non negligible probability of occurance. 

In particular, unless N is very large compared to Sg, the noise 
associated with a zero being received is different from the noise 
associated with a one. Therefore Fried and Schmeltzer's analysis 
can only be valid if the number of noise photoelectrons per pulse 
is much greater than the number of signal photoelectrons. In this 
case the two analyses should yield the same results. 

As a check on our results we have calculated the error prob- 
abilities assuming fourty noise photoelectrons per pulse, log 

amplitude variances of 0.0 and 0.05 and error probabilities from 
-1 -3 

10 to 3x10 . For these cases the computed power margin required 

to compensate for scintillation agreed with that reported by Fried 
and Schmeltzer to within one db. The residual error is due to the 
fact that even for fourth noise photoelectrons the condition of 
equation 7 is only approximately satisfied. 

On the basis of our analysis we conclude that Fried and Schmeltzer's 
results contain a fundamental inconsistancy ; namely, the number of 
signal photoelectrons required to obtain error probabilities less 

_3 

than about 10 is so large that equation 7 can not be satisfied 
unless the log amplitude variance is very small. Thus the assumption 
of background limited operation is invalid except for the case of 
weak scintillation or very high bit error rates. 
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OPTIMUM THRESHOLD 

The optimum threshold in the presence of scintillation is 
shown in figure 3 as a function of the average number of received 
signal photoelectrons per pulse. As is to be expected the optimum 
threshold decreases with increasing scintillation. For S Q greater 
than about 10 the curves are very nearly straight lines whose slopes 
are dependent on the strength of the scintillation. This linearity 
can be understood in the no scintillation case by noting that when 
eS is large compared to N equation 6 reduces to 

T = TrT(T7eT S ' (8) 

With scintillation we may replace equation 8 by a linear relation 
of the form 


T = a(C £ )S 0 + e (9) 

where the parameters a and e can be determined by linear least mean 
square fits to the data of figure 3. To a reasonable approximation 
we may take e to be 2.5, independent of C . The slope a, on the 

X/ 

other hand, is strongly dependent on as shown in figure A. The 
dependence of the slope on C can be represented by the relation 

Ai 

7.5X10 3 C 2 

a(C^) = . 1 45 ( e £ + 1) 


( 10 ) 
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Perhaps the most striking feature of figure 4 is the distinct 
knee that occurs for values of near 0.02, and corresponds to the 
transition from signal strength limited to scintillation limited 
conditions. For much larger than 0.02 the slope of the threshold 
curve varies only slowly with increasing C whereas for weak 

JO 

scintillation the slope is a very strong function of C^. In fact 
we can approximate a fairly well in these two regions by straight 
lines as shown by the dashed curves in figures 4. Hence 


V 


a i C £ S 0 + b i S 0 + C 


( 11 ) 


where the coefficients a^ and b. take on either one of two values 
depending upon whether is greater or less than the break point. 

A piecewise linear approximation of this sort is appealing from an 
engineering point of view since it provides a convenient model for 
implementing an adaptive threshold detection system. One could en- 
vision, for example, a detection system in which the mean and variance 
of the received signal power were continuously monitored and the 
threshold set in accordance with equation 11. Such a system would 
require only linear operations and a single discrete discontinuity 
to control the threshold level. 

As Titterton and Speck have pointed out (footnote 4 of reference 
2) this analysis is applicable to systems in which the threshold is 
varied on the basis of a long term estimate (ie. on the order of 
seconds) of scintillation and is limited to systems with bit rates 
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of 10 7 or less. Tycz, Fitzmaurice and Premo [5] have considered a 
perfectly adaptive threshold system for both the log normal and beta 
channels. 


SUB-OPTIMUM THRESHOLDS 

In practice the detection threshold of an optical communications 
system operating through the atmosphere will never be perfectly op- 
timized. If the threshold is fixed then changing atmospheric condi- 
tions will deoptimize the system and even if an adaptive threshold 
is used the system will be incapable of precisely tracking changes 
in signal strength and turbulence. In order to properly predict a 
communication system's performance it is necessary to know the ex- 
pected bit error probabilities for non-optimum detection thresholds. 
Investigation of the performance of non-optimum threshold detection 
is also necessary to the design of adaptive threshold systems. 

Figure 5 shows the bit error probability as a function of de- 
tection threshold with log amplitude variance as a parameter. The 
data plotted is for the case of 40 signal photoelectrons and no 
background photoelectrons per pulse and an extinction ratio of 15 db. 
Other choices of parameters yielded curves which were essentially 
similar. 

Inspection of figure 5 shows that with a fixed detection thresh- 
old a decrease in scintillation will always result in an improvement 
in system performance even though decreasing the scintillation de- 
optimizes the system. That is to say that the bit error probability 
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decreases more rapidly with decreasing log amplitude variance than it 
is increased by the corresponding deoptimization of the threshold 
level. Thus if one selects the detection threshold that is optimum 
for the strongest expected scintillation and predicts the error rate 
on this basis, one is assured that the system performance will not 
be worse under conditions of weaker scintillation. This strategy 
might be appropriate if one wishes to insure that a given error rate 
will be obtained under all expected conditions. 
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Figure 1. Probability of bit error for an on-off binary PCM 
optical communications system as a function of the number of 
signal photoelectrons for log amplitude variances between 0,0 
and 0.1. One background photoelectron per pulse and a modula- 
tion extinction ratio of 15 db. was assumed. 


Figure 2. The Loss Factor, defined as the additional trans- 
mitted power required to compensate for atmospheric scintil- 
lation, is given as a function of the number of received 
photoelectrons per pulse and the log-amplitude variance of 
scintillation. Other parameters are the same as the preceed- 
ing figure. 


Figure 3. The optimum detection threshold as a function of 
number of signal photoelectrons per pulse and log amplitude 
variance of scintillation. Other parameters are the same as 
in the preceeding figures. 


Figure 4. The slope of the linear portion of the detection 
threshold curves (figure 3) as a function of the log amplitude 
variance of scintillation. The dashed lines indicate a piece- 
wise linear approximation discussed in the text. 


Figure 5. Dependence of the bit error probability on detec- 
tion threshold for non-optimum detection. Error probability 
is plotted as a function of the normalized threshold, 1/S n , 
for S Q = 40. u 
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NEAR FIELD ANTENNA PATTERNS OF OBSTRUCTED 
CASSEGRAIN I AN TELESCOPES 

ABSTRACT 


The near field antenna pattern of an obstructed optical trans- 
mitter, such as a cassegrainian telescope, has been computed using the 
Fresnel approximation to the scalar diffraction integral. A number of 
cases have been considered. These include uniformly illuminated, ob- 
structed circular apertures, and apertures coherently illuminated with 
a divergent gaussian beam wave such as may be obtained from the TEM qq 
mode of a laser. Other beam profiles that may also been considered. 

The techniques used for the numerical evaluation of the Fresnel 
integral are discussed and results for the optical systems to be used 
on the Marshall Space Flight Center's High Altitude Aircraft Optical 
Communications Experiment are presented. 
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INTRODUCTION 


The Marshall Space Flight Center's High Altitude Aircraft Test 
for Visible Laser Optical Communi cations (AVLOC) will use cassegrainian 
telescopes as transmitting optics on both the air-to-ground and 
ground-to-air links, ^ For very long ranges, such as earth to sat- 
ellite, cassegrainian optics will cause no difficulty provided that 
collimated or focused beams are used. Under these conditions the 
shadow of the secondary mirror will be filled with diffracted light 
and the intensity maximum of the telescopes antenna pattern will lie 
on the optical axis of the telescope. That is to say that the far 
field diffraction pattern does not contain the shadow of the central 
obstruction. Many optical tracking and communication systems, including 
those employed on the AVLOC experiment, operate over much shorter path 
lengths so that both receivers are in the near field of their respective 
transmitters and hence see a doughnut shaped radiation pattern with a 
dark shadow on the optical axis. Even very long path lengths will not 
eliminate the central shadow if the transmitted beam has a divergence 
greater than the beam spreading due to diffraction by the aperture. In 
the case of divergent beams, one is in effect always in the near field. 

Since most optical communications systems require a divergent 
transmitted beam to keep the pointing errors less than the beam spread, 
it would appear that some technique must be used to eliminate the 
central shadow, or else cassegrainian optics (or any other obstructed 
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optical system) cannot be used. One might, for example, offset the 
transmitted beam with respect to the axis of the tracking system by an 
angle sufficient to remove the receiver from the dark central portion 
of the beam. This is undesirable, however, since the inclusion of such 
a point-ahead angle will complicate the pointing and tracking problems 
and is also wasteful of the transmitted energy. A more attractive 
solution would be to reduce the diameter of transmitted beam and dis- 
place it laterally from the optical axis of the telescope so that it 
uses a small transmitter aperture located in the unobstructed portion 
of the telescope. Other techniques, such as dithering the beam at a 
high rate to effectively smear it over a large area, or including 
diffusers in the optical system at appropriate points, might also be 
considered. 

In the case of the AVLOC test, there are two effects which could 
conceivably make the inclusion of a pointing offset or other technique 
unnecessary. The ranges to be employed in the experiment are such that 
the receiver is not located very deep in the Fresnel region. This is 
especially true of the downlink where the receiver will lie in the 
transition region between the near and far fields. At these distances 
there may be enough spreading of the beam by diffraction to provide 
usable levels of illumination near the axis. In addition atmospheric 
effects will tend to smear the beam and help to eliminate the central 
shadow. 

The primary purpose of this report is the compute the antenna 
patterns for the transmitting optics to be used in the AVLOC communi- 
cations system and to estimate the amount of illumination in the central 
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shadow of the cassegrainian obscuration* The theoretical relations 
and computer programs developed are quite general however, and can be 
directly applied to any obstructed optical system.. Although other 
investigators have studied the far field patterns of obstructed 
optical systems [2], to our knowledge the near field region has never 
been analyzed., 

We have computed the near field diffraction patterns for a 
cassegrainian telescope system transmitting a coherent gaussian 
beam* Numerical calculations have been made for a number of ranges 
and beam divergence angles typical of the conditions which will be 
encountered during the AVLOC test* Calculations have also been made 
for shorter ranges representative of the Marshall Space Flight Center 
Optical Range between Madkin and Bradford Mountains* These calcu- 
lations were performed for the purpose of comparison with beam profile 
measurements made by MSFC personnel on the ground transmitter system. 
In addition we have investigated the near field pattern of an 
obstructed telescope in which the beam profile has been redistributed 
by means of an optical axicon* The effect of atmospheric turbulence 
in helping to fill the central shadow is also considered* 



THEORETICAL BACKGROUND 


We may consider the diffraction pattern of the cassegrainian 
telescope as that of a circular aperture containing a central opaque 
disk and illuminated with a divergent gaussian beam.* Let U(x^,y^) 
be the complex representation of the field at an observation point 
(x q> y , z) and U(x^, y^) be the field at a point (x^ y^ o) within the 
aperture, z being the distance between the receiver and transmitter . 

We also define 

a - radius of the central obstruction 
b radius of the telescope objective 
A « wavelength 
k = the propagation number 
P = laser power 

R = radius of curvature of laser beam 
(0 = beam width 

0 = beam divergence (half angle) 

In the near field, the field distribution of a coherently illuminated 
aperture is given by the Fresnel approximation [3] 


The uplink beam is in fact not gaussian due to the inclusion of 
an energy redistribution device in the optical system. The effects of 
the beam redistribution will be considered separately „ 


4 



5 


[(x r x o )2 + <y r y o )2] | dx i dy i 

(i) 

where A f is the unobstructed portion of the aperture. We will assume 

that the laser is operated in the TEM mode so that U, (x^.y.,) is 

oo 1 V J 1 

given by the gaussian beam approximation 
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( 2 ) 


Ttbj 


0) 


222 2 1/2 
where , The factor (2P/iru> ) is included so that 

|u| will be the illuminance in watts per square meter. Substituting 

equation 2 into equation 1 and converting to polar coordinates we 

obtain 


U (r ,6 ) 
o o ’ o 


-aY 


2 n b 

/ / e xp [ jk 


2 -t 


2 

iru) o 


2R 


exp 



exp 


Jr + r ^ - 2r.r 005(0,-0 )] 

2z 1 o 1 o 1 o 


r l dr l d0 l 


(3) 


A phase factor -je has been dropped from equation 3 since it will 
not effect the intensity I = |uj . It should also be noted that 
contrary to common usage u is the distance at which the field, rather 
than the intensity, falls to 1/e of its maximum value and that the 
signs in the exponents have been chosen to make R positive for 
divergent beams instead of for focused beams. We now apply a well known 
integral representation for the Bessel Function, 
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J 0 (0 
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2tt 

J exp {-j 5 cos(0-<|i)} d0 
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to obtain 


u ( r ) 

o o 


2tt 

Xz 


VT / ex ^ { “r } exp{ ^ I ( i + i> r i 2} J o° 
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r i dr i 


Again a phase factor which will not influence the intensity has been 
omitted. 

Equation 5 is the Fresnel Diffraction pattern for an obscured 
aperture illuminated by a gaussian beam and is valid for all ranges 
likely to be of interest in optical communications. The persistence 
of the central shadow for large values of z can be understood by 
inspection of equation 5. Consider for a moment the case of a colli- 
mated beam. The radius of curvature of the beam is infinite and 
equation 5 reduces to 


U(r o ) 


2* /~2P 

Xz V 2 
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"I b 

/ 

a 
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exp{ — 5 —} exp{ j 
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2 z 
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r l dr l 


( 6 ) 


for z small U Q (r) must approach the geometrical optics limit of a 

gaussian beam with sharp central and peripheral shadows. However, 

jkr l 

for large z the term { — y — -} vanishes and equation 6 further reduces 
to the well known Fraunhoffer approximation 


U(r ) 
o 



exp{— y-} J r dr 

2 o z 11 


z 


( 7 ) 
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which has an intensity maximum at r = 0<, Thus the condition for the 

o 

elimination of the central shadow for a collimated beam is that 


z >> 



max a 


( 8 ) 


i.e^ , the range must be great enough to place the observer well into 
the far field „ Now for the case of an uncolliraated beam this con- 
dition becomes 


1 

R 



<< 



min. 


(9) 


If the beam is converging the Fraunhoffer condition is satisfied near 
the focus (R = -Z) • However, if R is small and positive, as it is for 
a divergent beam, the inequality of equation 9 cannot be satisfied for 
any value of Z, no matter how large* Hence one is in effect always in 
the near field and the central shadow will be observed at all ranges. 
For R small compared to Z equation 9 may be written as 


or 


R >> 



max 


R >> 



( 10 ) 


which corresponds to a beam divergence half angle of 


e < 

m 


_2_ 

kb 


( 11 ) 
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Here 0 is the largest beam divergence for which the central shadow 
m 

will disappear at long ranges. 

The preceding result may be interpreted physically as follows. 

For a divergent beam the geometric shadow of the central obscuration 
is a cone with it's apex located at the center of curvature of the 
wave front and a half angle 0. Light from the unobscured part of the 
aperture will be spread by diffraction; the angular divergence due to 
diffraction being of the order of X/b. If 9 is less than X/b then 
for large distances the light will eventually fill the central shadow, 
however, if 0 is greater than X/b the spreading due to diffraction 
is not great enough to compensate for the beam divergence and the 
shadow is never illuminated. 



ON AXIS ILLUMINANCE 


On the z axis equation 5 takes a particularly simple form since 
k r«. r 

here r = 0 and J ( ) = 1„ Then 

o z 


v°> 



r i dr i 


( 12 ) 


where 


a 



i „ k / 1 i 1 % 

+ 1 2 ( 5 + 


(13) 


to 

It will be useful to investigate the on axis illuminance since it 
will provide a check on the accuaracy of the numerical integration 
needed to evaluate the off axis Intensity. Equation 12 may be integrated 
directly to yield 


v°> • ! V 11 ? sr [ 

* 7TM y 


2 2 
+ab -aa n 

e - e J 


7T (0 / 

The intensity along the axis is then given by 


1(0) = 
1(0) - 


|u 0 (°) 


4z 2 1 a f 2 


2P 


7TW 


2 2 2 2 
/ ab aa s , a*b a*a v 
(e - e ) (e - e ) 


(14) 


(15) 

(16) 
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which reduces to 


1(0) k t ! . . ( e ( “ + “* )b2 + e ( “ + 

2ir z I a I w 


(aa^ + ot*b^) (a*a^ + ab^), 
— e — e J 


(17) 


'Rea(a^ + b^) /T ,2 , 2., , 

- 2e cos(Ima(a - b ))] 


(18) 


1 ( 0 ) 


2P 


tuo 2 [(!+§/ 


4z 

,2 *• 

k ti) 



e 





- 2e 


2 2 
a +b 


2 

cu 


cos 


r k J . l w 2 ,2. , 

[ 2 ( E + 7 )<a - b )] 


(19) 


The first exponential terms in equation 19 arise from the loss 
of the center portion of the beam from obscuration by the cassegrainian 
secondary and the loss of the tails of the beam at the objective. 

The harmonic term is due to the changing number of Fresnel zones con- 
tained within the aperture as one moves along the axis. This result 
is the alternating maxima and minima in the intensity along the axis 
characteristic of Fresnel dif fraction* 
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For the ranges to be used in the AVLOC experiment z/R »1 


and (2z/kto ) <<1. Equation 19 then simplifies to 


1(0, - -rfy 1 e- 2(b/ “ )2 + a- 2(a/ “ )2 

7TU) Z /R 

- -(a 2 +b 2 )/o, 2 ,k d lv .2 ,2.,, 

- 2e cos (y (-+ -) (a - b )}] 

( 20 ) 

But since 0 = o,/R (21) 


2 

and ir(z 0 ) is just the geometrical area of the beam in the receiver 
plane we may let I’ ( 0 ) denote the intensity which one would compute 
neglecting the central obstruction, the gaussian beam shape and all 
diffraction effects. Then 


1 ( 0 ) = I'( 0 ) 2 (e 2(a/ai) + e ~ 2 (b/w) 


„ -(a 2 +b 2 )/w 2 .k ;i . Is { 2 , 2 U 

_ 2e cos (y (r + 7 ) (a - b ) ) 


( 22 ) 


Clearly I (0) will have maxima and minima when the cosine term is ±1. 
Then 


*max/min * 


2 (e~ <a/uJ )± e - (b/ “ )2 ) 2 


(23) 


If we choose m = b and a = b/3 then I = 3.18 I' and I . = 0.55 I', 

max arm ’ 

or a factor of about 6 to 1 difference in intensity. If one is going 
to work on axis it would be worthwhile to choose the aircraft 
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range so as to take advantage of the maxima. The ability to do this 
depends, of course, upon the spacing of the maxima being great enough 
to allow the aircraft to remain near the maxima rather than fluctuating 
between maxima and minima. We may estimate this spacing by noting 
that maxima will occur when 


--- (| + |) (b 2 - a 2 ) = (2n + 1) ir 
n being an integer. Solving for z we obtain 


Z 

m 


(2n + 1) X 

,.2 2 , 

(b - a ) 



or 


AZ = 
m 


(2n+l)A 
(b -a ) 


~b 


(2n+3) X 
(b 2 -a 2 ) 




&Z - 


m 


Z - 
m 


[l/Z + : ■ > ■■■ 2 
m b 2 - a 2 


2X Z 


AZ = 
m 


m 


2X Z 


(b 2 - a 2 ) 


/ (1 + 


m 


,.2 2 ^ 

(b - a ) 


but 2X2/ (b 2 - a 2 ) << 1, therefore 
m 


2X Z 


2 


m 


AZ = 


m , 2 2 

b - a 


(24) 


(25) 


(26) 


(27) 

(28) 


( 29 ) 
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We choose a * 10 cm, b = 30 cm, X = 6328 A and Z = 20 km. AZ then 
is of the order of 6 km. so that the variation of I (0) with Z is 
clearly slow enough to allow the aircraft to remain near a maximum. 
Inspection of equation 25 shows that while the spacing between 
maxima is independent of R their position is not. Since R can be 
controlled by adjusting the beam divergence 6 it should be possible to 
locate the receiver near an axis illuminance maxima for both the 
uplink and downlink beams. 

The on axis illuminance which we have computed corresponds to the 
bright central spot usually observed in the Fresnel diffraction from 
an opaque disk [6]. The question remains whether or not this spot is 
broad enough to provide usable powers. In order to answer this 
question it is necessary to investigate the variation of illuminance 
with the distance off axis. 



ILLUMINANCE OFF AXIS 


The illuminance at a point off axis is given from equation 5 


2p v 2 2 2 

K V z) (s 2 + c 2 ) 

TTul) 


(30) 


where 


b -(r/w) kryr 2 

S = / e J ( -) • sin(0r- ) r.dr.. 

* O Z 1 11 

a 


(31) 


and 


b , , , 2 kr r „ 

C = / e j ( — — -) cos(6r. ) r-dr 

O Z 1 -L JL 

a 


(32) 


6 ■ I <5 + <33) 

The evaluation of S and C is relatively difficult since the argument 
of the Bessel function is such that neither the asymptotic nor ascending 
power expansions coverage rapidly. Furthermore, the trigonometric 
functions oscillate so rapidly that numerical integrating requires a 
very large number of points to accurately sample the function. This 
in turn requires that the integrand be evaluated with high precision 
to avoid an unacceptable accumulation of round off errors «, Several 
methods of evaluating C and S were considered, including an expansion 
by repeated integration by parts leading to a series of Bessel functions 
of increasing order and fixed argument* This series, which resembles 


14 



15 


the classical Lommel function expansion, was so weakly convergent that 
it proved to be even less efficient than direct numerical integration* 
The method finally chosen for the evaluation of S and C was a 
Simpson^ Rule integration with two correction terms* The zero order 
Bessel function was computed using an ascending power series [4] 


v« =E<- 1,r 


(f 0 2r 


<r ! ) “ 


(34) 


for values of £ less than 10 and Hankel's asymptotic expansion [5] 


J o (0 ~ {P (e> cos X + Q(£) sin x) ( 35 > 

£ ■> oo 


P 


1 - 


3 2 -5 2 


2! (80 2 4 ! (8C) 4 


(36) 


Q 



2 2 2 
1-3 »5 *7 

3»(8£) 3 


0 e 0 


(37) 


X - £ - tt/4 (38) 

for £ greater than 10.. Double precision arithmetic was used through- 
out. 

We have investigated the errors in the numerical integration, both 
by consideration of the error in the calculation of the integrand, by 
examining the magnitude of the correction terms in the Simpson f s Rule 
algorithm, and by . varying the number of points used in the integration. 

It is concluded that the maximum error will not exceed 0*1%. Furthermore, 
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it was found* that the value of 1(0) calculated from equation 30 agrees 
very closely with that found from equation 22 „ 

The computer programs used for these calculations are described 
in Appendix A» 



EFFECTS OF ENERGY REDISTRIBUTION 


In the previous sections the field distribution in both the ground- 
based and airborne transmitting apertures was assumed to be a truncated 
gaussian beam mode. The ground based telescope in fact will not have 
a gaussian intensity distribution in the exit pupil since it contains 
a device, called an axicon, which redistributes the beam in such a way 
that all the energy in the beam is incident upon the unobscured 
portion of the aperture,, For the purpose of estimating the order of 
magnitude of the illuminance near the axis the truncated gaussian 
approximation was adequate. It is desirable, however, to be able to 
compare the theoretical diffraction patterns with beam profiles 
observed over relatively short horizontal paths. In this way estimates 
of atmospheric spreading of the beam may be obtained. For this 
purpose it will be necessary to consider the effect of the energy 
redistribution of the near field pattern. 

Two possible physical realizations of the axicon are shown in 
figure 1. The axicon shown in Fig. la displaces the beam outward in 
a radial direction producing an energy distribution in the exit pupil 
usually referred to as a "displaced gaussian" distribution. The 
axicon illustrated in Fig. lb not only displaces the energy outward 
from the optical axis but also inverts the intensity distribution in 
the sense that light from the center of the entrance pupil is routed 
to the outside of the exit pupil and visa-versa. This type of 
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Entrance 

Pupil 


Exit 

Pupil 




a. plandromic gauss ian 



Entrance 

Pupil 


Exit 

Pupil 


b. displaced gaussian 


Figure 1. Two Axicon Devices 
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axicon produces an intensity dist rib ution" which has been described as 
"plalndromic gaussian". 

We assume that the energy redistribution device in the ground 
based telescope acts as a plalndromic axicon [7] i„e.. It maps the 
energy contained in an annulus of radius r* in the entrance pupil of 
the optical system into an annulus of radius r in the exit pupil as 
shown in Fig, 2a, The radii of the two annuli are related by 


i b ~ a i 

r = b = — r’ 


(39) 


and T is the truncation point of the original beam. This equation is 
equivalent to Peters and Ledger's equation (9) (reference 7) except 
that we have allowed an arbitrary truncation point for the original 
beam while Peters and Ledger have considered only beams which are 
truncated at a distance equal to the radius of the inner obscuration. 
Our results are therefore more general in that it will allow for an 
axicon with a lateral magnification different from unity. 

The redistributed beam profile is obtained from the mapping 
equation (equation 39) and from the conservation of energy. 


E (r)*27Trdr = E (r') "27rr*dr’ 
P g 


(40) 


where E^ and E^ are the plalndromic distribution in the exit pupil 
and gaussian distribution in the entrance, pupil respectively. Noting 


that E 


U we have 


iy r) 


(X^I) 

v rdr 


1/2 


|U g (r') 


(41) 
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in -u 1/2 T 2 o 

l U D (r) l ■ f b-r7 ,(_ i~ E) “ pr ' ( u (b - a)> (b ' r) ’ 


(42) 


Assuming that the beam wave fronts are approximately spherical 
we may add a phase term and a normalization constant to obtain 


1/2 1/2 T 2 2-, 

u (r) = C-^) 7—7 exp{ " — :v) (b ' r) } 

P TT CO 


w(b - a)' 


exp{j ~ r 2 } 


(43) 


Equation 43 is now substituted into equation 1 and following the same 
procedure as before we obtain the equivalent forms of equations 30-34, 
viz. 


2P , I , 2 ,k» 2 /. 2 ^ „ 2. 


«*.•»>■ V + V 

TT(i) 


(44) 


where 


b b - r -a 2 (b - r 1 ) 2 kr..r o 

S = / ( -) e J ( -) sin(Br 1 ) r-.dr. (45) 

t> J r, oz 1 J.X 

* a 1 


o z 


and 


. . 1/2 2 ,. .2 

b b - r -a (b - r ) k r 

C - } ( ) e T , 1 o 

p J r 1 J t “ 

r a 1 o z 


2 > 


J ( -) cos(6r 1 ) r t dr- (46) 

o z 1 11 


a and B given by 




a = — /(b - a) 

CO 


( 47 ) 
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Similar expressions can be obtained for a displaced gaussian. In 
this case equation 37 is replaced by the mapping 



(48) 


and the remainder of the derivation follows as before. The diffraction 


pattern from a displaced gaussian is then given by equation (44) with 


S and C replaced by and CL where 
p P J D D 


v ^ 1/2 2 . .2 , 

b a + r.. -o (a + r 1 ) kr 1 r _ 

S D = / E — ~ — ~] e J ( — ■—) sin(&r 1 ) (49) 

a 1 


and 


v ^ 1/2 2 , >2 

b a + r- -a (a + r.) kr.r 9 

C = / [ — - — -] e J ( — — 2-) cos (6r ) r.dr 1 (50) 

D ' r- 02 111 

a 1 


Numerical calculations have been made for a number of ranges and 
divergence using a plaindromic gaussian distribution and a displaced 
gaussian distribution. The numerical techniques employed were 
identical to those used for a gaussian distribution except for the 
choice of the input function. These results are presented in the 
following sections. 



RESULTS 


Numerical calculations have been made for the optical communi- 
cations systems which will be used by the Marshall Space Flight 
Center’s High Altitude Aircraft Optical Communications Experiment 
(AVLOC) • For the ground to air system the transmitter is a 24 inch 
(61 cm) diameter cassegrainian telescope with an 8 inch (20 cm) 
diameter central obstruction. The aircraft to ground system employs 
a 4 inch (10 cm) telescope with a 1 inch (2.54 cm) obstruction. 

The transmitted powers and wavelengths are 1 watt at 5880 nm. uplink and 
5 milliwatts at 632.8 nm, downlink. Ranges from 50,000 feet (15 Km) 
to 150,000 feet (45 Km) and beam divergence angles from 0 to 150 arc 
seconds were considered. The parameters for a typical. set of calcu- 
lations are shown in Table I. 

Table I 

Diameter (in) Beam Power Range Beam Divergence* 

System Aperture Secondary Watts 1000 ft. Arc Seconds 

uplink 24 8 1 80 0 

uplink 24 8 1 80 25 

uplink 24 8 1 80 50 

downlink 4 1 .005 80 25 

downlink 4 1 .005 80 50 

downlink 4 1 .005 80 75 


half angle 
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The assumed range (80,000 feet) corresponds to an aircraft 
approximately 10 miles from the ground station at an altitude of 
60,000 feet. For the purpose of these calculations we have assumed 
a gaussian beam whose beam width equals the outer diameter of the 
aperture. While this is not the optimum matching to the aperture, it 
is an adequate approximation, as the final results are not very 
sensitive to changes in the truncation point except for small changes 
in the overall amoliuudeo 


The results of these calculations are shown in Figs , 2-8 along 
with the geometrical optics approximation to the gaussian beam. The 
geometrical optics approximation was obtained by scaling the beam 
profile by the relation 


Z + R 

r = — r— r- 
o R j 


(51) 


For the uplink beam the diffraction pattern shows the charac- 
teristic monotonic decrease of intensity within the geometrical 
shadows, indicating that one will be deep within the Fresnel region 
as expected. Examination of the uplink beam patterns indicate that 
the intensities within the shadow may be expected to be quite low. 
Furthermore it seems unlikely that smearing by the atmosphere will 
be sufficient to produce usable illumination in the shadow. In each 
of these patterns the central maximum (Fresnel^-Arago bright Spot) is 
clearly visible. With the exception of the collimated beam it is 
far too narrow, however, to be of much use in a practical communi- 
cations system. 
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Figure 3. Beam Profile, Uplink. 
| Z - 80,000 feet 

I 0-0 arc sec. 
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For larger beam divergence angles , the solution tended to oscil- 
late very rapidly in the directly Illuminated region. This behavior 
was also observed at shorter ranges for smaller values of beam 
divergence. An example of this can be seen in Fig. 3 where the 
oscillation was so rapid that it could not be plotted. We have 
therefore indicated the limits of oscillation by the shaded area. This 
oscillation was at first a cause of considerable concern as we feared 
that it might be a spurious result due to error buildup in the com- 
puter calculations, A number of tests were made, however, which 
convinced us that the oscillation does correspond to real fluctuations 
of the Fresnel pattern. These tests included the following: 

1) Computing the intensity at very closely spaced points for 
one case. When this was done it was found chat the fluctuations were 
smooth oscillations rather than random variations which might be expected 
from error buildups. 

2) The number of points used in the Simpson's rule of integration 
was varied. It was found that the change in computed intensity due to 
changing the number of field points was much less than the magnitude 

of the oscillation. 

3) A gaussian integration algorithm was used to compute a few 
field points. These agreed closely with the results of the Simpson*s 
Rule integration. 

4) In all cases the magnitude of the oscillation was greatest for 
shorter distances of axis. This is the region where the integration 
algorithm should be the most accurate. 

5) No oscillation was observed in the shadow. 
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The rapid oscillation is of no importance as it will be obliviated 
by atmospheric effects a 

For the downlink the situation is somewhat improved* As can be 
seen from Figures 5-8 the intensity in the shadow oscillates, indi- 
cating that at the assumed range one is approaching the Fraunhoffer 
region. The central maximum has broadened and the intensity levels 
near the axis have increased. Since the calculated pattern will be 
further spread by atmospheric effects and by effects such as tracking 
system jitter it is expected that usable levels of illuminance will be 
found near the center of the pattern without additional beam spreading 
or pointing offsets • 

Results at ranges other than 80,000 feet were qualitatively 
similar* As would be expected the filling of the central decreased at 
shorter ranges and increased somewhat at longer ranges. . There was very 
little significant difference in the shape of the diffraction pattern 
at any of the distances considered however* Calculations made with 
plaindromic and displaced gaussian distributions displaced a somewhat 
different shape radiation pattern in the directly illuminated region 
but no significant differences in the fields within the shadows. 

We have also calculated the transmitter antenna patterns for the 
24 inch cassegrainian telescope at a range of 26,400 feet (8*15 Km), 
corresponding to the length of the Marshall Space Flight Center Optical 
Range between Madkin and Bradford Mountains. These calculations were 
made so that MSFC personnel could compare the actual antenna patterns 
of the ground based transmitter with the predicted pattern. Two repre- 
sentative patterns are shown in Figures 9 and 10. The beam divergences 
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are 12.5 and 25 arc seconds respectively and' a plaindromic gausslan 
beam was assumed in both cases. 
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Plalndromlc Distribution. 
Z » 26,400 feet 
6 * 25 arc sec. 
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atmospheric effects 


In a previous section it was suggested that smearing of the trans- 
mitted beam by atmospheric turbulence might be sufficient to eliminate 
the central shadow. In order to estimate the magnitude of the atmos- 
pheric effects we may adopt a ray optics approach and consider a 
differential" element of area in the transmitting aperture. A ray 
leaving this area will satisfy the ray equation 

d(ns)/d& = grad n (52) 

where n is the index of refraction of the atmosphere, s is a unit 
vector tangent to the ray, and di is a differential element of path 
length along the ray. Integrating equation 52 along the optical path 
we obtain 

_ z 

n's' - ns = / (grad n)d£ (53) 

o 

For large path lengths the deviation of the ray will not depend on the 
value of the index of refraction at the end points. We may therefore 
set n and n* equal to unity and obtain 

z 

s' - s = / (grad n)d£ (54) 

o 
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The angular deviation of the ray is just proportional to the vector 

difference s - s f . Following Chernov [9] we may obtain from equation 

53, after some algebraic' manipulations ,'the mean square deviation of 
2 

the ray LB . 


A 8 


= / / dr^d^ (grad^ grad 2 B(x 1 -x 2 ,y 1 -y 2 ,z 1 -z 2 


)) 


(55) 


B is the index of refraction correlation function which, for the 
Kolmogorov atmospheric model, is given by 

B(r) . j K ” r(8/3) slnW3) ^2 siair R -8/3 ^ (56) 

K 

o 

2 

where is the index of refraction structure constant and K q , are 
spatial frequencies corresponding to the outer and inner scales of 
turbulences (L q , -S, ) respectively. Equation 56 may be substituted 
into equation 55 and the integrals performed to give [10] 





(57) 


2 

Here and Jl have been assumed to be constants, as would be the 

case for a horizontal path. For vertical paths the appropriate 

variation of and with altitude would have to be included. 

We will use the horizontal path expression and assume the maximum 
2 

value of C^j and a minimum value of in order to establish an 
upper limit on the amount of scattering by atmospheric turbulence 
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which can be expected. Since I* o is much larger than L q we may drop 
the last term of equation 57. 

A0 2 » 5.7 C^ 2 Z Jl 0 " 1/3 (58) 

The energy from a single differential element will be spread 
over an area whose size is of the order of ajwhere a is obtained 
directly by taking the product of the mean square scattering angle and 
the distance from the scattering center to the observation point and 
integrating over the entire optical path. 

a 2 = f (z - S) 2 (^“) d£ (59) 

This integration yields 

a 2 - 1.9 ^ z 3 £ 0 “ 1/3 (60) 

Equation 60 expresses the variance of a rays lateral displacement 
in the plane of the observation point. We now wish to determine the 
probability that a ray originating at a given point in the transmitter 
aperture will pass through the observation point. Since the phase 
fluctuations of a wave passing through the turbulent atmosphere are 
known to be normally distributed, it is therefore reasonable to assume 
that the angular deviations of the ray is also normally distributed. 

This conclusion could also be reached by considering that for propagation 
distances that are large compared to the inner scale of turbulence, 
the net angular deviation of a ray is the sum of a large number of 
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small deviations o Hence' appeal' to the central limit' thereon leads 
to a normal distribution. v The lateral displacement of the ray is 
just the sum of the individual" angular deviations, weighted by the 
lever arm, i.e», weighted by the distance from the scattering point 
to the observation point. Therefore, the lateral displacement of 
the ray is a random variable which is the sum of a number of .inde-*- 
pendent gaussian random variables and hence should be normally dis- 
tributed. While the foregoing arguments may not be absolutely 
rigorous they do indicate the reasonableness of a normal distribution. 
Under any circumstances the error introduced by assuming a normal 
distribution should not be great. 

If we assume a gaussian distribution for the lateral ray dis- 
placements the contribution to the intensity at the observation point 
from a point (x^,y^) in the transmitting aperture is just. 


dl 


scatt 


lyviU 2 

(2tt) a 2 


2 

X 

exp { 2 } dx i d y;L 

2a 


(61) 


where x is the distance from the observation point to the point where 

the undeviated ray intersects the observation plane. If x Q y are the 

2 

coordinates of the observation point then x is given by 

x 2 = (x q - (1 + f) x ^ 2 + <y o " (i + f^) 2 (62) 


The total scattered intensity is then given by 


scat 


- // 


lyvi) 

(2tt) a 2 


exp {- [(x o " (1 + f )x l )2+(y o ~ (1 + f )y l> 
2 a 


2]}dx i dy i 


( 63 ) 
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where the integration extends over the unobseured portion of the 
aperture . 

To simplify our calculations let us assume an aperture ; 
illuminated by a beam of uniform intensity. We will also restrict 
ourselves to consideration of observation points on the optical axis 
of the transmitting telescope. With these restrictions equation 63 
becomes on changing to polar coordinates 


scat 


= -7 / exp { hr (1 


2a 


+ f > 2 


r l dr l 


(64) 


[ scat " (R/(z ' + R))2 [exp {-(1 + f) 2 > 

2 a 


exp {-(1+f-) 2 -hr 1 
R 2 a 


(65) 


Equation 65 may be evaluated numerically to give an estimate of 

the amount of filling of the central shadow. If we assume strong 

turbulence we may take C^ 2 to be on the order of 10 2 meters 2 and 

the inner scale of turbulence to be about 1 mm. With these values 

0 is found to be 1.5 m , If we then assume a 25 arc second beam and 

a range of 20 Km, we find that the intensity on the optical axis is 

about 2% of the intensity in the transmitting aperture. Thus under 

the assumed conditions the beam smearing would be appreciable. For 

-9 -1/3 

weak turbulence, however, might be as small as 10 meters , 
in which case the on axis intensity is only 


,_1_ -830 

^25 


) 



41 


times the intensity in the aperture. Clearly for weak turbulence 

the spreading of the - beam is entirely negligible. It should be 

remembered that these calculations apply only to a horizontal range 

since the altitude profile of C N 2 and SL q have been neglected. For 

2 

vertical paths the beam smearing will be much less since a 
decreases and i Increases with altitude. 



CONCLUSIONS 


The calculated beam profiles indicate that for the ranges and 
size optics to be used on the MSFC-AVLOC experiment the central 
shadow of the cassegrainian secondaries will be very dark and hence 
on axis operation will be impossible unless provisions are made 
to eliminate the shadow. Diffraction of light into the shadow 
was found to be negligble and the bright central spot, while 
present and strong, is far to small to be of any practical use. 
Smearing of the beam by atmospheric turbulence has also 
been considered but it was found that except for very long, horozontal 
paths and strong turbulence conditions the effect would be small. 
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APPENDIX A 


COMPUTER PROGRAMS 

1. Near Field Diffraction Pattern 

This program computes the near field diffraction pattern of a 
circular aperture or obstructed circular aperture by direct numerical 
evaluation of the Fresnel diffraction integrals A double precision 
Simpson’s Rule integration is used to evaluate the integrals The 
program as listed below will allow the aperture to be illuminated 
with a beam of arbitrary divergence angle and either uniform, gaussian 
or plaindromic gaussian. intensity profile „ A later modification will 
also accept a displaced gaussian beam. 

Input to the program consists of three data cards. The first 
card specifies the type of field distribution which is to be used, l.e., 
gaussian, plaindromic, or uniform, beginning in column one. The first 
character is compared with the characters ’G’ , 'P* and r U f stored in 
an array TCLASS and a branching flag is set and passed to the subroutine 
FRES*8 via the lab led common block /TYPE/. If the first character 
is other than 'G', ’P’ or *U* an error message is generated and 
execution terminates. The second data card contains the diameter of 
the obscuration (which may be zero) , the diameter of the aperture in 
inches, the beam power in watts and the wavelength in angstrom units. 

The card format is 4F10.4, The third data card contains the range, 
in feet, the beam width parameter in inches, the beam divergence angle 
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in arc seconds, the distance off axis of the first field point to be 
computed and the interval at which field points are to be spaced, 
both in feet, the number of field points to be , calculated and the 
number of points to be used in the Simpson's Rule integration. 

The input parameters are printed in the output listing and are 
then converted to metric units. The number of points specified for 
the Simpson's Rule algorithm is also checked and if it is even, it 
is increased by one. 

The subroutine DSMPC is called twice to compute the real and 
imaginary parts of the Fresnel Diffraction integral for the first 
field point. The intensity is then computed by taking the absolute 
square of the field. The off axis distance of the field point is then 
incremented and the intensity at the next field point computed. This 
process is continued until the specified number of field points have 
been calculated. 

After one profile has been calculated the program reads a new 
data card which may modify the range, beam width, beam divergence, 
field point and/or number of points in the integration. Execution 


continues until all data is exhausted. 
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// EXEC FORTGCLG,PARM,FORT='LIST,NAME=AVGO 
// COND.GOONLY SUPPRESS GO STEP 
//FORT.SYSIN DD * 

C PROGRAM TO COMPUTE FRESNEL DIFFRACTION PATTERN OF A CIRCULAR 
C APERTURE WITH A CENTRAL OBSCURATION 

C USING A DOUBLE PRECISION SIMPSON RULE WITH CORRECTOR TERMS 

C 

C 


c 

A 

s 

DIAMETER OF OBSCURATION 

INCHES 

c 

B 

ss 

DIAMETER OF APERTURE 

INCHES 

c 

LAMBDA 

= 

WAVELENGTH 

ANGSTROM UNITS 

c 

W 

B 

WIDTH OF GAUSSIAN BEAM 

INCHES 

c 

POWER 


LASER POWER 

WATTS 

c 

Z 

■ 

RANGE 

FEET 

c 

RHO 


DISTANCE OFF AXIS 

FEET 

c 

THETA 


BEAM DIVERGENCE 

ARC SECONDS 

c 

NPTS 

= 

NUMBER OF POINTS IN INTEGRATION 

L 

c 

T 

- 

TRUNCATION POINT FOR PLAINDROMIC FIELD 


C 

EXTERNAL FRES 

COMMON W,K,RHO,PI ,TWOPI ,R,Z,NTYP 
COMMON/TYPE/FASSA s PASSB ,t,wi ,IDIST 
DIMENSION LINE (80) 

INTEGER TCLASS(3)/'U’ ,*G' ,'P'/,TNUM/3/ 

REAL*8 A,B , POWER, LAMED A, PI ,TW0PI ,K,Z ,W , THETA, R,RH0 ,RH01,DRH0 ,RMAX , 
*SSM,SCM,E1»E2,INT 
*,PASSA,PASSB,T,WI 
PI = 3 d 141592653589 79 3D0 
TW0PI 2 =2o0D0 * PI 

C READ FIELD DISTRIBUTION TYPE AND SET BRANCHING FLAG 
READ(5 ,908) LINE 
IDIST=0 

DO 3 KK=1,TNUM 

IF(LINE(1) „EQ»TCLASS(KK)) IDIST=KK 
3 CONTINUE 

IF(IDIST.EQ.O) GO TO 10 
IF(IDISToNE „ 3) GO TO 11 
READ(5 ,911) T.WI 
WI=1„ 2700025 D-02 * WI 
T=l, 2700025D-02 * T 
11 CONTINUE 

C READ PARAMETERS AND CONVERT UNITS 
RE AD<5, 901) A,B, POWER, LAMBDA 
WRITE (6, 902) A, B, POWER, LAMBDA 
A=l<, 2700025D-02 * A 
B=l 0 2700025D-02 * B 
PASSA-B-A 



noon 
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PASSB=B 

K=2oODlO * PI / LAMBDA 

1 READ (5, 903, END-800) Z »W, THETA, RH01, DRHO ,NRHO,NPTS 
WRITE (6 ,904) Z,W, THETA 
WRITE (6, 909) LINE 
W-l. 2700025D-02 * W 
Z=0s 3048006096D0 * Z 
THETA = 4o 84814D-06 * THETA 
R = W/THETA 

RHOl - 0. 3048006096DO * RHOl 
DRHO = 0. 3048006096DO * DRHO 

PARAMETERS HAVE BEEN CONVERTED TO METRIC UNITS 
NOW CHECK NUMBER OF POINTS *** IT MUST BE ODD 

N-MOD (NPTS ,2) 

IF(NoEQ.O) NPTS-NPTS+1 
WRITE (6,906) NPTS 
DO 2 J-l.NRHO 
RHO - RHOl + (J-l)*DRHO 
NTYP-1 

CALL DSMPC(SSM,FRES,E1,NPTS,A,B) 

NTYP-2 

CALL DSMPC(SCM,FRES,E2,NPTS,A,B,) 

C COMPUTE INTENSITY 

INT= ( (K/Z) **2) * ( (2 . 0D0*P0WER) / (PI* (W**2) ) * (SCM**2+SSM**2) 

C PRINT RESULTS 

2 WRITE(6,907) RH0,INT ,E1 S E2 
GO TO 1 

10 WRITE (6, 910) 

800 STOP 

901 FORMAT (4F10 » 4) 

902 FORMAT ('I* , 2 5X, 'PARAMETER LIST’ //10X,'A= * ,F10. 5 , 1 INCHES ' ,5X, ’B-' , 

*F10.2,* INCHES' /10X, 'POWER-’ ,F6„2,' WATTS * /10X, 'WAVELENGTH- 1 ,F8<,2, 
** ANGSTROMS ’ ) 

903 FORMAT(5Dl4o6,2l5) 

904 FORMAT (10X, 1 RANGE-* ,F12.2,* FEET’ /10X, ’ BEAM WIDTH-’ ,F10. 2, 1 INCHES 
*’ /10X, ’BEAM DIVERGENCE-’ ,F10.2, t ARC SECONDS') 

906 FORMAT (10X, 'USING A ’,15,’ POINT SIMPSONS RULE ’) 

907 FORMAT ('O’ ,15X,’ DISTANCE OFF AXIS «’ ,F10. 4, /5X,' INTENSITY =’, 

*D24„ 16/25X, 'ERROR TERMS' ,5X,2D26. 16) 

908 FORMAT (80A1) 

909 FORMAT (10X, 'TYPE DISTRIBUTION = ’,80A1) 

910 FORMAT ( ’ ILLEGAL FIELD DISTRIBUTION SPECIFIED EXECUTION TERMIN 
*ATED ' ) 

911 FORMAT (2F10„ 4) 

END 
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The subroutine DSCMP is a double precision Simpson's Rule Algorithm 
with two correction terms . The correction terms are printed to provide 
an estimate of the accuracy of the integration routine. 


SUBROUTINE DSMPC(SUM,FUN,E ,N,FLL,FUL) 

DOUBLE PRECISION SUM, SI, S2»S3,S4,FN,FLL,FUL,DX X,FJ,Y,E,FUN 

51 = 00.0D0 

52 = 00.0D0 

53 = OOcODO 

54 = 00.0D0 
FN = N-l 

DX = (FUL-FLL) /FN 
DO 10 1=1, N 
J=I-1 
FJ=J 

X=DX*FJ +FLL 
Y=FUN (X) 

IF(IcEQ„l»OR„IcEQ,N) GO TO 1 
IF(I.EQo2.0R.IoEQ 0 (N-l)) GO TO 2 
K=M0D(I ,2) 

IF(K) 3,4,3 

1 S1=S1+Y 
GO TO 10 

2 S2=S2+Y 
GO TO 10 

3 S3=S3+Y 
GO TO 10 

4 S4=S4+Y 
10 CONTINUE 

SUM=(S1+4„0D0* (S2+S4) + 2c0D0*S3) * (DX/3.0D0) 

E=(-4. ODO^Sl + 7„0D0*S2 + 8. 0D0*(S4-S3) + FUN (FLL-DX ) + FUN( 

1 FUL+DX)) * (DX/(-9oODl)) 

SUM=SUM+E 

RETURN 

END 

/ /LKED c SYSLMOD DD DSNAME=UA0L02.EEG00118.PWRFLW(AVLOC2) ,DISP=SHR, 

/ / UNIT=DISK , V0L=SER=UA0L02 , SPACE= 

//GO.SYSIN DD * PROGRAM CHECKOUT DATA 

GAUSSIAN 

8.0000 24.0000 1.0000 4880.0000 

80.000000D 03 24.000000D 00 1.500000D 02 1.000000 00 1.000000D 00 1 31 
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The subroutine FRES*8 is called. by. DSMPe to compute the integrand. 
Two Bessel function approximations are used; a Jacobi asymptotic 
expansion and an ascending power series^ The program selects the 
appropriate routine each time depending’ on the value of the argument. 


REAL FUNCTION FRES*8(R) 

REAL* 8 W,K,R 9 RHO,Z ,ARG,A,A2 ,A3 ,A4,A5 »P ,Q,PI ,TW0PI,RND ,ARGJ, JZRO, 
*TRIG,RAD 
*,A6,A7,FAC,XI 

* , U , ARGE , ARGRT ,P AS SA ,P AS SB , T 
*,WI 

COMMON W,K,RHO,PI,TW)PI,RAD,Z,NTYP 
COMMON /TYPE /PAS S A , PAS S B ,T , WI ,ITYP 
C COMPUTE ZERO ORDER BESSEL FUNCTION FOR LARGE ARGUMENT USING 
C JACOBI ASYMPTOTIC EXPANSION 

C 

ARGJ=K*R*RHO/Z 

IF(ARGJ.LT.l.ODOl) GO TO 1 

A=8. 0D0*ARGJ 

A2=A*A 

A3=A*A2 

A4=A*A3 

A5=A*A4 

A6=A*A5 

A7=A*A6 

P=1.0D0 - 4.5DO/A2 + 4.59 375D2/A4 r- 1.50077812 D05/A6 

Q=1.0D0/A - 3, 75D1/A3 + 7. 441875D3/A5 - 3.62 3071875D06/A7 

RND-ARGJ/TWOFI 

NR=RND 

RND=NR 

RND=ARGJ-RND*TWOPI 
RND » RND-PI/4.0D0 

JZR0=DSORT(2.0D0/(PI*ARGJ))*(P*DC0S(RND) + Q*D IN (RND)) 

3 CONTINUE 

ARG=0 . 5D0*K* ( 1 » 0D0 / RAD + 1.0D0/Z)*(R**2) 

IF(NTYP.EQol) TRIG=DSIN (ARG) 

IF(NTYP»EQ„2) TRIG=DCOS (ARG) 

GO TO (11,12,13) ,ITYP 

11 FRES=TRIG*JZRO*R/l. 414213562373095D0 
RETURN 

12 CONTINUE 

FRES - DEXP (— 1 . 0D0* ( (R/W) **2) ) * TRIG * JZRO 8 
RETURN 

13 ARGE=— 1. 0D0* ( ( (PAS SB-R) /PASSA)* (T/WI) ) 

ARGRT= (PASSB-R) /R 
IF(ARGRT.LE.O.ODO) GO TO 14 

FRES=DEXP (ARGE) *JZRO*R*DSQRT (ARGRT) *TRIG*W*T/ ( I*PASSA) 

RETURN 
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14 FRES=0 . ODO 
RETURN 

1 A= -l.ODO * (ARGJ / 2 . ODO) **2 
AL=1 . ODO 

FAC=1 . ODO 
JZR0=1. ODO 
DO 2 1=1,20 
XI=I 

FAO=FAC*XI 

A1=A1*A 

A3=A1/(FAC**2) 

A4=DABS (A3) 

IF (A4 . LT . 1 . ODO-16 ) GO TO 3 
JZRO=JZRO + A3 

2 CONTINUE 
GO TO 3 
END 


2. Program for Calculation, Field Intensity on Axis 

This program evaluates the Fresnel diffraction integral for 
field points on the optical axis. In addition, it computes the 
field distribution from the geometrical optics approximation. 


PI-3, 14159 

4 READ (5,1) A, B,W,P, WAVE, THETA, Z 
WRITE (6, 2) A, B,W,P, WAVE, THETA, Z 
A=0.0127*A 
B=0 . 0127*B 
W=0. 0127*W 
WNUM= 2 . OE 10 *P I /WAVE 
Z=0. 304 8*Z 

R=W/(4. 848E-06*THETA) 

C=2.0*P/(PI*(W**2)) 

D= (1. 0+Z/R) **2 + ((2.0*Z)/(WNUM*(W**2)))**2 

C=C/D 

X=(A/W)**2 

Y=(B/W)**2 

ARG= (WNUM/2 . 0) * (A*A-B*B) * (1 . 0/Z+l. 0/R) 

PD=C* (EXP (-2 . 0*X) + EXP (-2. 0*Y) - 2.0*EXP(-1.0 (X+Y))*C0S(ARG)) 
WRITE (6, 3) PD 
WRITE (6, 13) 

RAT=(R+Z)/R 

XW=W*RAT 

C0= (2 . 0*P) /PI* (XW**2) ) 
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XLp-A*RAT 

XHI=B*RAT 

WRITE (6, 11) XLO , XHI 
DX= (XHI-XLO) / 10 . 0 
DO 10 1=1,11 
X=(I-1)*DX + XLO 
PD=CO*EXP (-2 . 0* (X/XW) **2) 

10 WRITE (6,12) X,PD 

11 FORMAT (10X, 'LIMITS OF GEOMETRICAL SHADOW ’,25 4.6) 

12 FORMAT (2 OX,' X= 'E14.6,' INTENSITY ' ,E14.6) 

13 FORMAT (’0 GEOMETRICAL OPTICS APPROCIMATION' ) 

GO TO 4 

1 FORMAT (6F10. 4, F10.0) 

2 F0RMAT('0* ,10X,6F10.4/11X,F10.0) 

3 FORMAT (10X, 'INTENSITY = ',E12.4,' WATTS PER METER SQ ') 
END 
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STUDY OF ATMOSPHERIC EFFECTS ON LASER COMMUNICATIONS SYSTEMS 


ABSTRACT 


The atmospheric scintillation has been observed at 10.6pm over a 
3.2 Km path using a heterodyne optical communications system. The 
observed probability density functions agreed reasonably well with a 
log-normal model for small receiving apertures but showed noticable 
departures from log-normalicy for apertures above 4 cm. The power 
spectral density of the scintillation has also been computed. 
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CHAPTER I 


INTRODUCTION 

The atmosphere is an Inhomogeneous non-isotropic media which is 
usually in a state of turbulence. The atmosphere is characterized by 
its temperature, wind velocity, and humidity. The variations of these 
parameters comprise non-stationary random processes and as a result the 
fluctuation of the index of refraction of the media Is also a non- 
stationary random process. Electromagnetic radiation at optical wave 
lengths propagating through the atmosphere will be greatly affected by 
the random fluctuation in the index of refraction. This causes a very 
complicated scattering to occur which results in amplitude and phase 
variations of the wave. Clearly, these fluctuations are also of a random 
nature. 

The distortion induced by the atmosphere on the propagating wave is 
of great concern in the development of optical tracking and communication; 
systems since performance can be seriously degraded due to this effect. 
For example, amplitude fluctuations cause the signal-to-noise ratio to 
be reduced in . incoherent detection systems. Atmosphere distortion tends 
to reduce the coherence of the wave which in turn reduces the effective 
power level of the signal. This will also decrease the signal-to-noise 
ratio. Loss of coherence is a problem in systems utilizing coherent de- 
tection where heterodyne action must be achieved. Wave front tilt due 
to phase variations of the wave induces errors in optical tracking re- 
ceivers since they view this tilt as an apparent angle of arrival. 

An important aspect in the design of optical systems that is de- 
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pendent on a knowledge of atmospheric distortion is the size of the optical 
components. Atmospheric effects can be somewhat overcome by employing 
large receiving optics which will tend to average out the wave fluctuations. 
This advantage is limited by the high expense and difficulty in fabrication 
of such components. The designer must then choose optimum size components 
which require that he have a thorough knowledge of the atmospheric pro- 
blem. It is clear that there is a pressing need for an accurate 
mathematical model of the atmosphere. 

Statistical methods must be employed to analyze this problem since 
the wave fluctuations are random processes. The first step in developing 
a statistical model is to determine the probability density function of 
the random process. Theoretical considerations have predicted a log- 
normal distribution for the amplitude and normal for the phase. This 
theory has been experimentally verified for visible wave lengths, but re- 
sults of current investigations in- the infrared region of 10.6 microns 
have been inconsistent. 

D. L. Fried^ has made scintillation measurements at this wave length 
over a 1 km. range using a point detector. His results do not confirm 
the hypothesis that intensity scintillation is log-normally distributed. 

He suggests that this may be a genuine feature of 10.6 micron scintilla- 
tion but draws no definite conclusion since detector noise and nonlinearity 

problems in taking measurements could have influenced his results. 

2 

Richard Kerr of the Oregon Graduate Research Center has conducted 
multiwave length laser propagation studies over a mile path and claims 

O 

confirmation of log-normal statistics for wave lengths of 4880A and 10.6 

3 

microns. In addition Fitzmaurice, Bufton, and Minott have also con- 
cluded that scintillation at 10.6 micron fits the log-normal model. Their 
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work was done over a 2.4 km. path. Both investigators used point 
source detection. 

The effect of the atmosphere on 10.6 micron propagation is important 
since the popular CC^ laser emits radiation at this wave length. This 
type laser is attractive for application in optical systems due to its 
high efficiency and high output power capability. In addition, the 
atmospheric effect at this wave length is much less than that at visible 
wave lengths. 

The purpose of this study is experimentally to investigate the 
statisitical properties of scintillation and the signal-to-noise ratio 
of heterodyne detection for a CC >2 laser beam propagated over a 3.2 km. 
path. Both scintillation and heterodyne measurements have been made for 
a variety of receiving aperture sizes ranging from two to ten cm. 

A brief discussion of the theory which is referred to in current 
literature is presented in Chapter II. The necessary statistical con- 
cepts are introduced before a qualitative description of atmospheric 
turbulence is given. Finally, the physical significance of aperture 
averaging is discussed. 

Chapter III gives a detailed description of the experiment. De- 
scribed is the equipment, its alignment and check out as well as a 
discussion on the techniques used to make the measurements. 

The handling and reducing of the data is given in Chapter IV. This 
includes a discussion on the conversion of analog data to digital form 
for direct use on a digital computer. An outline of the computer pro- 
gram which reduces the data is presented. The theory used to calculate 
aperture effects is also given. 

Chapter V is concerned with interpreting the reduced data to de- 
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termine if the hypothesis of the log-normal distribution for intensity 
scintillation is valid for this wave length. This chapter also in- 
cludes results of calculations for the refractive index structure con- 
stant with and without aperture averaging corrections. 

Chapter VI contains the summary and conclusions of this study as 
well as recommendations for further study. A complete documentation 
for the computer program is given in the Appendix. 



CHAPTER II 


THEORY 


A. Statistical Concepts 

It is necessary to give a discussion on pertinent statistical con- 
cepts as a prelude to presenting a qualitative discussion on the 
theoretical aspects of the atmospheric problem. 

The random processes are described in terms of parameters which are 
random variables. The value of any such function at a fixed instant of 
time is a random variable having definite probability density function. 
The process may further be described by its auto-covariance function at 
times and t ^ 

AC{f(t 1 ), f(t 2 )} = <ff(tp - <^f(tj> ][f(t 2 ) - <^<t 2 5>l> 2-1 

where indicates an ensemble average. The auto>-covariance function 

reduces to the correlation function 

BtfCtp, f(t 2 )] = <f( tl )f(t: 2 J> 2-2 

for processes where the mean value is zero. The auto-^covariance function 
characterizes the mutual relation between the fluctuations at different 
instants of time. The mean value of the random variable can be a con- 
stant or can change with time. Similarly, the auto-covariance function 
can either depend only on the difference between the times t^ and t ^ or 
else it can depend on the positions of the points on the time axis. The 
first case would occur when the statistical relation between the 
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fluctuations of the variable at different instants of time does not 
change with time, A random function is called stationary if its mean 
value doe 3 not depend on time and its auto-covariance function depends 
only on the difference between observation times. 

The mean value of the meteorological parameters of the atmosphere 
such as temperature, wind velocity, and humidity undergo comparatively 
slow and smooth changes. These variables are non-stationary processes 
if the definition of stationarity is strictly applied. It is difficult 
to determine which changes in the fluctuation are to be regarded as slow 
changes in the mean and which are to be regarded as slow fluctuations of 
the function. 

To avoid this difficulty and to describe random functions which 

have the above characteristics, the structure function is used instead 

of the correlation function. This function was first introduced by 
4 5 

Kolmogorov ’ , The basic idea behind this method is to use the difference 
function 


F (t) = f(t+i) - f(t) 2-3 

T 

instead of the non-stationary function f(t). For values of t which are 

not too large, slow changes in the function f(t) do not affect the 

value of the difference function which means that it can be considered a 

stationary random function. The function f(t) is called a random function 

with stationary increments. To derive an expression for the structure 

function consider the transformation of the correlation function for 

F(tJ and F(0: 

T i T * 


B<'l. * 2 > ■ < F T< t l>V t 2>!> 


2-4 



7 


B(t^, t 2 ) = ^f(tj+T) - f (t^) ] [f (t 2 +x) - f(t 2 )f> 


Using the algebraic identity 

(a-b) (c-d) = | [(a-d) 2 + (b-c) 2 - (a-c) 2 - (b-d) 2 ] 2-6 


we have 



Thus the correlation function is expressed as a linear combination of 
functions of the form 


w 



< 




- £<t J )] 
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which is called the structure function of the random process. The form 
of the structure function more commonly used 


D f (t) = 




X 


[f(t+T) - f <t) 3 



2-9 


is the basic characteristic of a random process with stationary in- 
crements. The value of D(t) characterizes the intensity of those 
fluctuations of f(t) which are smaller than or are comparable with t. 


B. Nature of Atmospheric Turbulence 

The statistical theory of turbulence was initiated in the works 
of Friedmann and Keller 6 . This theory was greatly amended in 1941 
when A. N. Kolmogorov and A. M. Obukhov 6 established the laws which 
characterize the basic properties of the microstructure of turbulent 
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flow at very large Reynolds numbers. The following discussion is drawn 
from Tatarski's^ work on the Kolmogorov theory. 

Consider the atmosphere to be a viscous fluid in a state of 
laminar flow. This flow can be characterized by its viscosity v, velocity 
v, and the characteristic length L. The quantity L characterizes the 
dimensions of the flow as a whole and arises from the boundary con- 
ditions of the fluid dynamics problem. This laminar flow will be stable 
if the Reynolds number 


v 

does not exceed a certain critical value. 

Suppose that for some reason a velocity fluctuation occurs in a 
region of size t of the initial laminar flow. The value of Reynolds 
number will increase and the laminar motion will lose stability. The 
result of this instability is the formation of a secondary flow or eddies 
within L which will have their own Reynolds number R^. As the Reynolds 
number for the overall flow is increased, will increase causing the 
secondary flow to break up into smaller scale eddies. These new eddies 
now give energy to even smaller eddies and the process continues until 
an eddy with a Reynolds number less than the critical value is formed. 

The atmospheric turbulence can be considered as consisting of many 
circulating eddies having different flow characteristics. Eddies are 
usually described in terms of an inner and outer scale of turbulence. 

These are measures of the characteristic sizes of the smallest and largest 
eddies which exist at the time of interest. Figure 1 may aid in visualizing 
this process. The outer scale is the physical dimension of the largest 
eddy. The inner scale of turbulence is roughly the size of the smallest 
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stable eddy. It can be more precisely defined in terms of a characteristic 
of the longitudinal velocity structure function 

D rr = < \ tV(7 l ) " V(7 2 )] ^ 2-11 

where V(rp is the projection of the velocity at the point along the 
direction of r, and V(r^) is the same quantity at the point 

2-12 
2-13 

For r on the order of L the structure function saturates. The inner 

scale of turbulence £, is then defined mathematically as the value for 

D where the functions in equations 2-12 and 2-13 intersect, 
rr 

Each eddy or cell in the field of turbulence can be considered 
locally isotropic and homogeneous, and as a result it will have a certain 
index of refraction, which we will assume to be constant throughout the 
cell. The index of refraction will in general differ from cell to cell. 

As an electromagnetic wave passes through each cell two things occur: 

First, the phase of the wave is advanced or retarded in a random manner 
due to the index of refraction of the cell. Secondly, the wave is 
scattered due to the interfaces between the cells. This causes the in- 
tensity of the beam to be distributed at random after traveling through 
a large number of cells. 

From the model of the atmosphere just developed some insight as to 
the statistical characteristics of phase and amplitude variation can be 


r = r 4* r 
2 1 


For r << l 


D - ar 
rr 


and for r >> SL 


_ 2 2/3 

D = c r 
rr 
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gained. Since the random variations in phase add to each other, the 

Central Limit Theorem^ can be applied to predict that the distribution 

of the phase across the wave front is normal. The amplitude of the wave 

has a distinctly different character. * After each refraction the intensity 

is the product of the intensity before scattering and the variation due 

to refraction. The intensity variations are then due to a product of 

probabilities. If the Central Limit Theorem is applied to the logarithm 

of the intensities, they will be normally distributed. This leads to a 

log-normal probability density function for the amplitude distribution. 

Because of this, it is customary to describe waves propagating through 

the atmosphere in terms of their phase cf»(r,t) and log-amplitude L^(r,t), 

where L (r,t) is given by 
a 




£n { 


A(r,t) | 


2-14 


or in terms on intensities 




2-15 


where A and I are mean values. Using these equations the complex repre- 
sentations of the wave becomes 


A(r,t) exp [L(r , t) + j<t>(r,t)J 


2-16 


C. Aperture Averaging 

Collection of light with large aperture optical systems tends to 
average out atmospherically induced intensity and phase fluctuations. 
This causes a smaller variance for both and alters the statistical prop- 
erties of the intensity variation. 
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Before discussing the principles of aperture averaging it is 
necessary to define correlation distance Consider the intensity at 

two points separated by a distance of r. For r equal zero the correlation 
function will be unity. As r increases, the correlation function decreases 
with zero as its lower bound. How rapidly the correlation function de- 
creases with increasing r, is related to the strength of atmospheric 

turbulence, r is defined as the distance at which the intensity variations 
o 

of the two points are no longer correlated. r^ could also be defined as 

the distance at which the variations become statistically independent. 

r usually varies from a few millimeters for very strong turbulence to 

several centimeters for mild turbulence. 

Consider the aperture shown in Figure 2 to be the aperture of an 

optical system which collects and focuses light from a diverging beam. 

Let the light intensity across the illuminated aperture be divided into 

n finite circles of radius r . The intensity within these circles will 

o 

be assumed to be highly correlated. The collection and focusing process 
can be thought of as adding the intensity contribution of each circle 
on the aperture in the focal plane. Averaging of the intensity fluctua- 
tions of each circle across the aperture will occur at the focus if the 
diameter of the aperture is much larger than such that it contains 
many circles. 

The intensity fluctuation in the circles of radius r^ are log- 
normally distributed. Since the aperture adds the intensities, it will 
also add the probability density functions. The Central Limit Theorem 
can then be applied to predict that the intensity variations at the focus 
should be normally distributed. In order to apply the Central Limit 
Theorem, we must assume that the average intensities of the circles are 
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of the same order of magnitude across the aperture and that the variance 
from circle to circle should not change significantly. 



CHAPTER III 


EXPERIMENTAL 

A. Description of Equipment 

The experimental measurements were made during the Summer of 1969 
at the Marshall Space Flight Center’s optical range located on Redstone 
Arsenal near Huntsville, Alabama. 

The transmitting system was located in an astronomical observatory 
on the crest of Madkin Mountain. The receiving and recording systems 
were located in the Astrionics Laboratory complex in a special building 
equipped with a large mirror periscope so that the optical equipment 
could be conveniently placed on the ground level yet have a clear optical 
path to the mountain, 'ihe height of the periscope was about 15 feet 
above ground level. The optical path extended in a southwesterly 
direction for a distance of 3.2 km. The transmitter was about 220 meters 
above the receiver so that the optical path was at an angle of 4° with 
the horizontal. Except for a few small buildings and a parking lot 
paved with bituminous material, the optical path lay mostly over wooded 
terrain. 

The transmitter and receiver were constructed by Minneapolis 
Honeywell Corporation for Marshall Space Flight Center and have been 
described in the literature . The transmitter consisted of a 5 watt C0 2 
flowing gas laser with a 10 cm. off axis cassegrainian collimator as 
shown in Figure 3. The laser was designed to have good short and long 
term frequency stability. This was accomplished in part by constructing 
the cavity of the low expansion material cervit, which has an expansion 


15 






17 


coefficient of less than 1 x 10 ^ (1/°C). To further stabilize the 
laser, water held at a constant temperature to better than 0.1°C was 
continuously circulated through the cervit yoke. The laser in the trans- 
mitter was frequency modulated by applying ^he modulation voltage to a 
piezoelectric cylinder on which one of the end mirrors was mounted. The 
other end mirror consisted of an Irtran output coupler that was also 
attached to a piezoelectric cylinder, which provided laser transition 
selection, A DC bias was applied to the cylinder to select the desired 
transition. In addition, the transmitter included a mechanical chopper 
that was originally intended to be used for alignment purposes. The 
transmitter unit also contained the necessary electronics to produce a 
modulation voltage for both carrier or direct modulation. A block diagram 
of the transmitting unit is given in Figure 4. 

The receiving unit was housed in a cabinet identical to that of the 
transmitter as shown in Figure 5. The receiver consisted of a 10 cm. 
off axis cassegrainian telescope, a local oscillator laser Identical to 
that of the transmitter, combining optics, and a mercury doped cadmium 
telluride detector. This is an alloy detector having a spectral re- 
sponse between 8 and 14 microns. The detector was cryogenic and required 
an operating temperature of 77°K, which was obtained by using liquid 
nitrogen. The operation of the receiving unit can be described with 
reference to Figure 6. The transmitter and the local oscillator signal 
are made spatially colinear by means of a germanium beam splitter and 
combined on the surface of the mercury doped cadmium telluride detector. 
The local oscillator frequency is offset by 10 Mhz from the transmitting 
laser. The 10 Mhz beat frequency produced by the detector is amplified 
with a 10 Mhz center frequency, intermediate frequency amplifier. This 
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21 


amplifier has a 2 Mhz bandwidth and a 110 db gain. The intermediate 
frequency amplifier is followed by a limiter that eliminates amplitude 
variations. A 10 Mhz discriminator provides an error signal for the 
local laser feedback loop which consists of a low pass filter and a 
feedback amplifier. The feedback amplifier drives a piezoelectric cylinder 
on which one end mirror of the laser is mounted. This provides auto~ 
matic frequency control of the local oscillator laser. 

The data acquisition system was located at the receiver terminal 
and consisted of an Ampex 14 channel analog F.M. tape recorder, a variable 
gain AC amplifier with good low frequency response, and two oscilloscopes 
used for monitoring purposes. Figure 7 gives a block diagram of this 
system. A spectrum analyzer was also available to check the 10 Mhz beat 
signal in the receiver. 

The monitoring of both the input to the amplifier and the recorder 
was necessary to insure that they were operating within their linear 
range. Especially critical was the input level to the recorder, since 
its linear range for input voltage was ±1 volt. To be safe we operated 
within ±.5 volt. 

B. System Alignment 

It was necessary to measure the laser and amplifier noise of the 
system to insure that it would not have a significant effect on 
measurements made through the atmosphere. Noise measurements were made 
with the transmitter and receiver placed a few feet apart with the local 
laser in the receiver turned off, so that only noise contributions from 
the transmitter laser, detector and receiver electronics would be 
present. The system was then operated in the same manner over a 100 
meter path through an enclosed tunnel one meter in diameter. The noise 




Figure 7. Block Diagram of Data Acquisition System. 
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of the local oscillator laser was determined by operating it into the 
detector in the absence of an incoming signal. The noise in all cases 
was found to be sufficiently low so that it would be negligible compared 
to the expected variation due to atmospheric scintillation. The linearity 
of the mercury cadmium telluride detector was determined by noting 
changes in DC voltage output for different power levels. The laser 
power was measured on one side of the germanium beam splitter using a 
Coherent Radiation Laboratories power meter and the DC variation in the 
detector was measured by a digital voltmeter. For incident power level® 
less than 300 mw. , the detector output was found to be linear. 

Alignment of the system over the 3.2 km. path proved to be a 
difficult task. Our first attempt was to bore sight a 60 power tele- 
scope mounted on the transmitter case, with the output beam. The idea 
was to aim the transmitting unit on the mountain at the laboratory 
periscope. This method worked over the 100 meter tunnel quite well, but 
the bore sight became misaligned when the transmitter was transported to 
the mountain. The second method employed to align the system involved 
the use of two visible lasers. The transmitter unit, now mounted on the 
observatory telescope stand on Madkin Mountain, was aligned by placing 
an argon laser directly in place of the receiving unit in the laboratory. 
The bright beam could easily be detected by the eye at the transmitting 
terminal. The position of the argon laser was adjusted until its beam 
was intercepted by the objective of the transmitter unit. The optical 
system of the transmitter was then adjusted so that the visible laser 
was focused onto the output aperture at the transmitter laser. Using 
heat sensitive paper as a position indicator for the infrared beam, the 
visible light and the invisible beam were made to coincide at two points 
in the optical system. Alignment of the receiving unit was accomplished 
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in a similar manner except that a visible laser could not be mounted in 
the same position as the now aligned transmitter unit. A small helium 
neon laser was mounted with a telescope on a tripod and located as near 
to the transmitter as possible. The telescope was then bore sighted to 
the helium neon laser. The laser-telescope arrangement was pointed by 
locating the top periscope mirror at the laboratory. A comer re- 
flector located at this mirror enabled a more precise aim as the re- 
flection of the red light could be seen with the telescope. The 
visible beam was then focused onto the detector in the receiver. With 
minor adjustments to the transmitter laser mount, close alignment was 
attained for the system. 

C. Measurement Procedure 

Scintillation measurements were made with the receiver laser 
inoperative so that only light from the transmitter and from the sun’s 
reflection off the observatory dome were intercepted by the receiver. 

This background light was cause for great concern since accurate 
measurements of the variations of the laser light could not be made in 
its presence. Since we could not filter out this unwanted light, it was 
necessary to record it. To accomplish this scintillation measurements 
were made by chopping the transmitter beam at 90 Hz by means of the 
mechanical chopper located in the transmitter. Figure 8 shows a sample 
of this chopped signal. 

Scintillation measurements were made in 90 second segments. A written 
log was kept on each data run giving the date and time of day, analog tape 
number, the channel number, aperture size, and weather conditions such 
as temperature, wind velocity and humidity. The intensity signal for 
each data run was read on the analog tape into marked time slots. One 
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Figure 8. Waveform of Signal Recorded for Scintillation Measurements* 
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of the channels was marked every 90 seconds with an audio tone code so 
that approximately 33 slots were available on each data channel. When 
sufficient signal was being received, data runs were made in succession 
using five different size apertures. The aperture size varied from two 
to ten cm. in diameter. The runs were made in succession to insure that 
atmospheric conditions remained constant over each set of five runs. 

Since the signal level for small aperture sizes was below that for 
the larger apertures, it was necessary to adjust the AC amplifier from 
one run to another using the monitor oscilloscope. Also, it was 
important to keep the dewar on the detector filled with liquid nitrogen 
so that the detector temperature would remain constant. 

Measurements of the signal to noise ratio and distortion were made 
on the communication system operating over the optical range. The signal- 
to-noise ratio of the heterodyne action was measured at the receiver by 
extracting the 10 Mhz beat note between the received signal and the local 
oscillator after it had passed through one stage of amplification. The 
signal was detected with a simple diode circuit and the resulting 
voltage recorded by the data acquisition system. These measurements were 
also made for different aperture sizes. 



CHAPTER IV 


DATA REDUCTION 

A. Analog-to-Digital Conversion 

All measurements were recorded on analog tape. In order to process 
this data it was necessary to convert it to a digital form suitable for 
computer input. The very large quantity of data collected necessitated 
electronic conversion to digital magnetic tape which could then be read 
directly into the computer. 

The first step in the conversion from analog- to-digital form con- 
sisted of recording a timing signal on a reserved channel of the analog 
tape. Eight standard time signals are broadcast by Marshall Space Flight 
Center’s Computation Laboratory on a frequency of 226.5 Mhz. The time 
signals are subcarrier multiplexed with pilot frequencies between 2.3 KHz 
and 70.0 KHz. The signal chosen was a rectangular pulse with a one- 
second repetition rate broadcast on a subcarrier frequency of 52.5 KHz. 
This signal is designated as 100/1000 AMR-D-5, and is coded with the 
time of day in Greenwich Mean Time by pulse width modulation. This timing 
signal had no relation to the actual time the data was recorded. 

The analog tapes were then read into a cathode-ray-oscillograph. 

Four data channels, the timing channel and the marker channel were re- 
corded simultaneously. The oscillograms were inspected and the sections 
of the tape to be digitized were selected. At this time bad data was 
identified and eliminated. The starting and stopping time for each time 
interval selected was read from the timing channel and recorded on the 
Computation Laboratories instruction forms . The analog-to-digital con- 
verter system was set to convert only those time intervals selected by 
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reading the timing channel. Approximately 30 one-minute segments were 
selected from each channel. A9 the running time for the analog tape 
was about 45 minutes , about two-thirds of the total data recorded on a 
tape channel was converted. 

The analog-to-digital conversion was performed by Marshall Space 
Flight Center's Computation Laboratory using an Astrodata type con- 
verter. Five data channels were digitized simultaneously. The five 
channels were sampled alternately at a rate of 5000 samples per second 
which amounted to sampling each channel at 1000 samples per second. The 
resulting binary coding was recorded on seven channel digital tape in a 
multiplexed format. 

The digital tape format consisted of a ten bit binary word so that 
2048 levels were available to represent analog signal levels between -.5 
volts and +.5 volts. Since seven track tape was used, each sample re- 
quired two tape characters consisting of a ten bit word plus the sign 
bit, blank bit and two parity bits. The data was recorded in records of 
2004 tape characters. The first four tape characters contained the time 
as read from the timing channel at which the first sample was taken. The 
remaining 2000 characters contained 1000 sample points, 200 from each 
channel, alternating between channels. Subsequent records were written 
until the segment was completed, at which time an end of file mark was 
written; therefore, each file on the digital tape contained one time 
slice or data run from the analog tape. In addition, the first record 
of each file was an identification record of 24 tape characters which con- 
tained the tape number and other information. 

B. Description of Computer Program 

A program was written for an IBM-360-50 computer to reduce the data 
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stored on magnetic tapes. This program reads data from the tape, changes 
the binary format to a fortran compatible form, then calls its various 
subroutines to perform the analysis. The principle problems which were 
encountered in writing this program concerned formatting the data for 
the computer and extracting the actual light intensity signal from the 
modulated square wave. When the signal was extracted it was either 
stored for spectral analysis or is used to construct a histogram. A 
Fourier transform subroutine or statistical subroutine was then called 
to analyze the signal. 

Development of a routine to extract the desired signal proved to 
be somewhat difficult since the sampling rate during digitization could 
not be accurately synchronized with the period of the square wave. The 
sampling rate of 1 KHz and the chopping rate of 90 Hz should yield 
approximately 10 samples per cycle of the square wave. In actuality, 
the number of samples per cycle varied between ten and twelve due to the 
sampling rate not being an integral multiple of the square wave frequency. 
The routine was designed to determine whether a particular data sample 
was a base point (i.e., from the part of the square wave when the laser 
beam was blocked by the chopper) or a signal point (when power was being 
received from the laser beam) » The problem was further complicated by 
the fact that the rise and fall times of the square wave were non- 
negligible so that about one percent of the data points were sampled 
during the switching transient and should be neglected. In addition, 
some of the data contained an occasional noise spike which should be 
eliminated. It was decided that the elimination of these spikes would 
not adversely effect the validity of the analysis so provisions for 
eliminating them were also included in the program. 
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The extraction routine operates basically as follows! A pre- 
processing routine reads the data from the magnetic tape and stores a 
record containing 200 data points into a common array. To begin the 
analysis twenty data points from the array are selected and their maximum 
and minimum value computed. Two limits, and are then set by the 
relations 
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where A and A . are the maximum and minimum values of the first 
max min 

twenty points, and P^ and are constants between zero and one half. 
Since the signal was inverted when it was recorded on the analog tape, 
the base line is greater than the signal, hence a particular point 
greater than is considered a base point, if it is less than it is 
considered a signal point. Points lying between and are assumed 
to be from the transient portion of the wave form and are neglected. The 
routine takes each point successively and determines if it is a base 
point, a signal point or neither. As a preliminary to processing, the 
first twenty points are scanned and the beginning of a base line segment 
of the wave form is found. Then new limits are set on the next 15 points 
and they are scanned and grouped into three arrays, a base line segment, 
a signal segment and a second base line segment. Each array may contain 
up to ten points. At this time another routine is used to compute the 
signal amplitude of the square wave for the group of signal points (as 
will be described later). The second group of base points is transferred 
into the first array, new limits are set using the next ten data points, 
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and a new group of signal points and base points are found to fill the 
second and third arrays, and finally their amplitudes computed. This 
process is continued until the 200 points from the first record have been 
used. At this time the routine pauses while the next record is read in. 
Processing then continues until the number of records called for (up to 
300) have been processed. 

This routine also contains several checks to handle possible 
irregularities in the data. During the search for either base or signal 
points if more than ten consecutive points are found, the routine will 
request that the next record be read in, which means the remainder of 
the bad record is discarded. Also if the number of unsuccessful scans 
while in the base or signal search phase exceeds ten, the routine will 
enter an error recycle phase. In this phase the routine skips 20 points 
and resumes processing as if it were at the beginning of a new run. If 
the error recycle phase is entered five times in a given record, a re- 
quest for the next record will be executed. When a new record is re- 
quested due to an irregularity in the record being processed the routine 
again treats it as it does the first record of a new run. In addition 
to the above, if three or less base or signal points are found in a 
given search, the error recycle phase is entered. 

During processing, a record is kept of each irregularity encountered 
and this information is printed in tabular form when the processing of 
a run is completed. If an excessive number of irregularities occurs in 
a given run, the results of that run must be suspected. 

The three arrays containing base and signal points found by the 
extraction routine are passed into a routine that determines the actual 
amplitude of the signal. Three methods for computing the amplitude were 
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tried: The first method took to base line for a group of signal points 

as the average of all the points in the group of base points preceding 
it and the ones following it. This is equivalent to considering the 
background light during the time when the laser beam intensity was being 
recorded to be the average of the background light recorded during the 
half-cycle immediately preceding the signal and the half-cycle immediately 
following it. The difference between the signal point and the average 
base line was taken as the amplitude of the laser beam at that instance. 

The second method of amplitude calculation considered was to re- 
construct the base line by fitting a least-mean-square curve to the base 
points. This method produced erratic results and also required additional 
computer time and was abandoned. 

In the third approach, the difference in the first signal point in 
a group and the last base point preceding it is taken as the signal 
amplitude. The difference in the last signal point in the group and the 
first base point following it gives a second amplitude. This method 
yields only two amplitudes per cycle but has the advantage that they are 
evenly spaced. Since the difference in time between signal and base 
points is very small, this method eliminates the problem of the unknown 
base line over the signal interval. 

The computer program contained both the first and third methods of 
signal amplitudes calculations. The method to be used was selected by 
a parameter read during execution. The values computed by this routine 
were either stored in an array to be used in the spectral analysis, or 
used to construct a histogram. Histograms are sometimes referred to as 
being the probability density function for discrete variables. 

The final segment of the program consists of the analysis routines. 
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The statistical routines accept the histogram for the intensity 
fluctuation which has been generated and computes the scintillation 
statistics. The routine computes and lists the corresponding value of 
the log-amplitude as defined by 

L ± = i £n [I ± /Y] 4-3 

where L and 1^ are the log -amplitude and the intensity for the ith class 
interval and I is the mean intensity. The frequency for each class and 
the cumulative probability are also listed. The routine also calculates 
the mean, standard deviation, skewness, and kurtosis for both the in- 
tensity and log-amplitude distribution. A chi-square test that checks 
the intensity distribution for a normal fit and the log-amplitude for a 
log-normal fit is also included. Appendix A includes a typical computer 
print out of the statistical analysis. 

The program includes routines to perform a spectral analysis on the 
intensity scintillation. In calculating the scintillation spectrum 

Q 

2 (8192) values of the beam intensity were extracted from the raw data 

using the preprocessing routines described previously. These values 
represent the intensity fluctuations of the beam sampled at 180 Hz rate. 
The spectrum of the sampled data was computed using the Cooley-Tukey 

9 

algorithm for fast Fourier transforms. The resulting spectra were dis- 
played by means of a plot routine. In cases where data points are dis- 
carded due to irregularities in the extraction routine, the omitted points 
were replaced with zero values. The purpose of this was to preserve the 
phase relationship of the signal, which is important in the integration 
operation of the Fourier transform. 
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C. Program Check and Parameter Adjustment 

To test the program several records were read from a magnetic tape 
and printed out for inspection. Each sample was classified either as a 
signal point, or a base point, or as a point from the transient portion 
of the waveform. This classification was purely subjective, yet in 
inspecting the data there was usually no questions as to how a point 

should be categorized. The same data was then read into the computer. 

* ~ 

A print out of all the pertinent variables at each decision branch of the 
extraction routine was obtained. The program was then run several times 
varying the parameters P^ and in equations 4-1 and 4-2 between runs 
and the results compared with the subjective analysis. Data having 
irregularities was also processed and the results compared to determine 
whether or not a segment should be omitted. Using these comparisons as 
a criteria, the parameters and P^ were set at 0.05 and 0.10 respectively. 
Therefore a point within 5% of the maximum base point or 10% of the 
minimum signal point would be retained while points between these limits 
were discarded. 

The statistical routines were checked out with sets of numbers which 
had a known probability distribution function. A routine readily available 
in the IBM Scientific Subroutine Package for the IBM 360 computer was used 
to generate a large set of normally distributed numbers. These numbers 
were then processed by the statistics routines in the same manner as the 
intensity fluctions were to be processed. The part of the routine which 
tests for the normal distribution fit gave very positive results, and 
that for the log-normal fit gave negative results. The test routine was 
then altered to generate numbers with a log-normal distribution. The 
-r*>«ults of the statistical analysis of this data strongly indicated a 
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log-normal fit. A description of the log-normal generator is given 
in Appendix B. 

D. Atmospheric Structure Constant and Aperture Averaging 

An important parameter in the statistical model for atmospheric 

studies used in current literature is the refractive index structure 
2 

constant C^, defined as a constant of proportionality in the relation 

D n = < \ tn l< r l>- n 2 (r 2 )l / > " C n ^ 

r = r 2 ~ r ± 4-4 

where D is called the structure function and is a measure of the 
n 

deviation of the index of refraction at two points separated by a 
distance r. is actually a measure of the strength of the turbulence. 

Fried^ gives a relation involving for an infinite spherical wave 

propagating a distance z in a turbulent atmosphere. The relation is 

C.(0) = .124 K ?/6 Z ll/6 c 2 4-5 

SL n 

where C. (0) is the log-amplitude variance. Equation 4-5 can be used 

2 

directly to obtain by noting that the standard deviation of the 
log-amplitude distribution which we have computed is the square root of 
0 ,( 0 ). 

A finite receiving aperture has the effect of averaging the in- 
tensity fluctuating from various parts of the wave front thereby re- 
ducing the variance of the scintillation. 

The effect of aperture averaging can be allowed for by using re- 

11 12 

lations developed by Fried 9 viz. 
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a s 2 = [f D 2 ] 2 6 CjCO) 4-6 

Where a g is the signal variance which corresponds to the square of the * 
standard deviation of the intensity fluctuation, D is the diameter of 
the receiving aperture, and 8 is an aperture averaging factor given by: 
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H(p/D) is the optical transfer function of a circular aperture 
H(p./D) = cos 1 (p/D) - p/D [1-Cp/D) 2 ] 1 ^ 2 4-8 


and C^(p) is log-amplitude covariance given by: 
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In the last expression a and b are the expansion coefficients 

n n 

for the modified confluent hypergeometric function and are given by 
the recursion relations 


a = -a - {(2n - 23/6) (2n - 17/6)/ (2n-l) (2n) } 
n n— I 

b = -b . { (2n - 11/6) (2n - 17/6) (2n-l)/(2n) (2n+l) 2 } 
n n-1 


4-10 
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where 

a = 1 b = 6.84209 

o o 

The intensity variance 0^.(0) can be related to the log-variance by: 

C T (0) = I q 2 [exp(4C £ (0)) - 1] 4-11 

2 2 

Equation 4-11 specifies C (0) in terms of a and I . I and a are 

1/ s o o s 

the mean and standard deviation of the recorded intensity. Combining the 
above equations we have: 





ttD [exp [40^(0) ] - 1] 
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where f(kp/4z) is the summation given in equation 4-9. Since a g /I^ is 
an experimentally determined constant, equation 4-12 is an integral 
equation for (0) . A special computer program was written to solve 
equation 4-12, The technique used is to evaluate the integral in 

equation 4-12 for a number of trial values of (0) using a fourth order 

a 

g 

Runga-Kutta integration. This gives a table of — as a function of 

o 

C (0) . From this table the value of C (0) corresponding to the measured 
o 

g 

— is determined using Lagrange-Hermite interpolation formula, 
o 


value of 



CHAPTER V 


RESULTS 

A. Probability Density Function for Intensity Scintillation 

It has been customary in the literature to test the hypothesis of 
log-normality of scintillation data by plotting the cumulative prob- 
ability function of the log-amplitude against a "probability scale" such 
that if the data is log-normal the resulting curve will be a straight 
line. The same method can be applied when testing the intensity amplitude 
for a normal distribution. Since this test would require considerable 
time if it were applied to every run, it was necessary to use a test 
that could be incorporated into a computer program in an efficient manner 
in order to determine which distribution each run more closely fit. 

The necessary statistical parameters to make this test were calcu- 
lated as described in Chapter III and were available on punched data 
cards. The skewness coefficient was chosen as the parameter to indicate 
the type of distribution. The skewness coefficient will ideally have 
zero value for perfectly normal or log-normally distributed data; however, 
for real data we expected a small value but somewhat greater than zero. 

We have chosen the skewness coefficient as the measure of closeness 
of fit in preference to the chi-square test since the chi-square depends 
upon the number of class intervals in the sample while the skewness 
does not. To use the chi-square on a comparative basis would require the 
generation of a table of chi-square for all possible number of class in- 
tervals and would lend to undesirable complexities in the computer pro- 
gram. 
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In order to use the skewness as a criteria for categorizing the 
data it was necessary to set bounds on its value* Since the lower 
bound is zero, it was only necessary to determine the largest value that 
the skewness could attain which would represent a suitable fit* This 
was accomplished by plotting several graphs of the cumulative probability 
for runs with values of skewness ranging from .02 to .79* The chi- 
square test made on each run was examined to insure that the overall 
results of the analysis were consistent* A selected number of these graphs 
are given in Figures 9-17* An inspection of these plots will show that 
the cumulative probability curve becomes nonlinear with increasing 
values of skewness. An inspection of many such plots indicated that the 
curve deviated from linearity much faster when values of skewness be- 
came greater than 0.15* For values from 0.0 to 0.15 the curve remained 
approximately linear * 

The results of the statistical analysis of each run was categorized 
in the following manner: If the skewness for the log-amplitude data is 

less than or equal to 0.15, the run was categorized as being log- 
normally distributed* If the value of the skewness meets the above re- 
quirement for the intensity data the run was considered to be normally 
distributed. If both skewness coefficients were greater than 0.15, the 
run was put in a "neither” category. In addition we have had to in- 
clude a category for those distributions which were sufficiently close 
to both a normal and a log-normal distribution that we could not 
distinguish between them. These runs which have approximately the same 
skewness for both distributions have been designated as "both". 

The results of the categorization for 196 runs are shown in Table 
I. Only runs with a small number of preprocessing irregularities have 
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APERTURE 

SIZE 

(CM) 

NUMBER OF RUNS 

TOTAL 

LOG 

NORMAL 

NORMAL 

NEITHER 

BOTH 

2 

4 

2 

2 

0 

0 

4 

17 

13 

2 

1 

I 

6 

30 

10 

6 

13 

s 

8 

43 

| 

1 2 

14 

! 1 

6 

10 

n 

33 

2 3 

37 

9 


Table 1. Results of Statistical Analysis of Scintillation Measurements. 
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been tabulated. From this table we see that the two cm. aperture runs 
were divided equally between normal and log-normal distributions. The 
results for this aperture size are not considered significant since there 
were only four runs with small irregularities, and these were taken with 
a very weak signal. The four cm. runs, which had a reasonably small 
aperture, yet the recorded signal was strong, show a strong tendency 
toward log w normality . As the aperture increases the number of runs which 
differ from log-normal also increases. This behavior may be due to the 
effect of aperture averaging. Since this process is additive in nature, 
it would cause the distribution to tend toward a normal curve according 
to the Central Limit Theorem. This combined with the difficulty in 
distinguishing normal data from log-normal data could lead to results of 
the type exhibited by the larger aperture runs. 

B. Calculation of Atmospheric Structure Constant 

The refractive index structure constant was calculated using 
equation 4-5 and taking the measured value of the log-amplitude variance 

C £ (0) . This equation does not allow for aperture averaging effects. 

—8 1/3 —81/3 

Values of C obtained ranged between 1.6 x 10 M and 8.7 x 10 M 
n 

for a group of runs where the aperture was varied rapidly from two to 
ten cm. The structure constant computed from two cm. aperture data was 
usually larger than that computed from ten cm. aperture data by a factor 
of from two to four. The time required to collect the data for such a 
group of runs was less than ten minutes. As a comparison, several 
groups of runs were made in a similar time period holding the aperture 
constant. The observed variation in the structure constant for these 
runs was usually about ten percent. This clearly indicates that the 
observed variation in the structure constant was due to aperture 
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averaging effects and not to changes in the strength of turbulence be- 
tween runs. 

Using the techniques described in Chapter III, section D, calcula- 
tions for the structure constant have been refined, allowing for the 
effects of aperture averaging. For the 212 segments of data which were 

analyzed, the values of the corrected structure constant ranged from 
-7 1/3 “7 1/3 

5.8 x 10 M to 9.0 x 10 M . These values are characteristic of 
very strong atmospheric turbulence, which agree with our subjective 
observations of the scintillation while the data was being taken. 

Table II shows four typical sets of runs for various aperture 
sizes. As can be seen, the variation in the structure constant is 
significantly reduced when the effects of aperture averaging are in- 
cluded. 

It should be noted that the inclusion of aperture averaging effects 

has a tendency to overcorrect the structure constant variation. This 

could be an indication of a systematic error in recording the data. 

Another possibility is that the deep scintillation conditions under which 

the experimental data was collected produced saturation effects which have 

not been considered. On the other hand the validity of the basic approach 

to the problem of aperture averaging in terms of the structure function 

13 

has been questioned . Therefore it is possible that the aperture 
correction we have used is not valid. In any case, this method of com- 
pensating for aperture averaging effects is a good approximation since 
the structure constant is far more consistent when corrected than when 
aperture averaging effects are neglected. 

C. Scintillation Frequency Spectrum 

The frequency spectrum has been computed using the Fast Fourier 
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Transform techniques described above. This calculation was performed 
only on selected data runs due to the time required for such a calcula- 
tion, Also the computed spectra were very consistent so it was felt 
that analysis of additional runs would yield little additional informa- 
tion. 

The computed spectra cover a range of frequencies from DC to 90 Hz, 
the upper limit being set by the 180 sample per second sampling rate. 

A typical spectrum is shown in Figure 18. Although the computer plotted 
frequency components out to 90 Hz only components out to 20 Hz are shown 
in the figure for the sake of clarity. Above 20 Hz the spectrum con- 
tinued to decrease linearly so that no appreciable frequency components 
above 40 Hz were observed. 

D. Heterodyne Detection 

The effect of atmospheric turbulence on the performance of the 
equipment operating as a heterodyne communication system was investigated. 
This was accomplished by recording the amplitude of the 10 Mnz 
heterodyne signal at the output of the I. F. amplifier. Also the trans- 
mitter was modulated with a 1 Khz signal which was recorded at the out- 
put of the receiver T s F. M. descriminator . It was found that neither 
signal showed any effect attributable to atmospheric turbulence large 
enough to be accurately measured. Even under conditions of deep 
scintillation encountered during the course of this experiment, the 
atmospherically induced noise was of the same magnitude or smaller than 
system noise. It was also found that clear voice communications were 
possible over this range under the worst conditions of scintillation 
encountered. 

While these results clearly indicate the feasibility of using a 
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heterodyne CC^ laser system to communicate through the earth atmosphere, 
they were somewhat disappointing since they did not permit a quantitative 
measurement of noise induced by atmospheric turbulence. 



CHAPTER VI 


SUMMARY AND CONCLUSION 

The results of the scintillation measurements made on the 10.6 
micron wave length laser beam tend to confirm the log-normal model for 
small receiver apertures* The data for the larger apertures did not 
seem to fit the log-normal or normal models with any consistency. One 
possible explanation for this could be aperture averaging. It is 
possible that the larger aperture sizes were not large enough with re- 
spect to the correlation distance to cause the distribution to be normal, 
but yet large enough to cause the distribution to deviate from log- 
normally. This in addition to the difficulty in distinguishing between 
the two distributions could have caused the results to be inconclusive. 

The value of the refractive index structure constant computed was 
found to lie within the range of values for this constant as calculated 
by Fried. The value of this constant was found to decrease with in- 
creasing aperture size. Equations developed by Fried were used to correct 
the structure constant for aperture effects. This technique seemed to 
give a slightly larger value of the aperture averaging effect than was 
observed. Although this may indicate an inaccuracy in the theoretical 
expression, a systematic error in the experimental measurements cannot 
be firmly ruled out. These calculations were significant in that they 
indicated in a quantitative manner the nature of aperture averaging. 

The spectral analysis indicated that low frequency components of 
the scintillation were predominant. The magnitude of the scintillation 
decreases linearly with increasing frequency. Above 20 Hz the 
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scintillation is negligible. 

The feasibility of optical heterodyne communication at 10 microns 
through extreme turbulence was demonstrated. To our knowledge, this was 
the first system of this kind to be operated through the atmosphere at 
this path length. We feel that the successful operation of the 
communications system under extreme scintillation conditions was a 
significant result. 

It would be a great advantage in this type experiment to reduce the 
data immediately after it is taken. Information gained from the speedy 
reduction should give the experimenter a knowledge of how the system is 
performing and aid in making better measurements. 

As a result of performing this experiment and surveying the results 
of other investigations, it is clear that the atmospheric problem is far 
from being completely solved. The log-normal model needs to be further 
verified for other wave lengths. The variance and the structure constant 
should be investigated under as wide a variety of weather conditions as 
possible and should be correlated to the variations of the meterological 
parameters. 

After having carried out this type experiment, the need for several 
refinements in the procedure was realized. Before any data is taken, a 
thorough analysis of system noise should be made. Sensitive noise 
measuring instruments should be employed. The noise of the transmitting 
laser should be recorded simultaneously with scintillation measurements. 
The background light effect should be further studied and perhaps a 
method other than signal chopping used to handle the problem. Mechanical 
stability of equipment is a problem that needs investigation since even 
very small vibrations could cause the beam to shift out of alignment 
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with the receiver. 

Aperture averaging effects need to be investigated with many 
different size apertures for all wave lengths. The relationship between 
correlation distance and aperture size needs to be determined. The 
determination of correlation distance in itself would be an interesting 


experiment . 
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APPENDIX A 


Introduction 

This appendix further describes the computer program employed to 
reduce the atmospheric data* The program is written in Fortran IV 
language for the IBM-360 Model 50 computer at the University of Alabama. 

The program operates in the following manner: The MAIN or supervisory 

routine accepts instructional and operating data for the program. It 
reads the atmospheric data from the magnetic tape and calls subroutine 
EXT to extract the intensity signal. EXT calls on subroutines LIMIT, 

FILL, HIST, and AMPX to perform the extration. When a file of data has 
been processed MAIN calls subroutine PRINT to print a table of the 
irregularities that occurred during extraction. MAIN then calls sub- 
routine STAT and/or FFT to perform the statistical and spectral analysis. 
STAT calls on subroutine CHI to perform the chi-square test which in turn 
calls on a Simpsons rule integration subroutine SIMP. FFT calls on two 
package subroutines FOURT and PLOT which perform and plot the spectral 
analysis . 

Routine Main 

Routine MAIN’s first step is to read the identification record from 
the magnetic tape. This is done by calling subroutine RID. RID 
actually does the reading and stores the data in a common array to be 
printed by MAIN in the next step. The numerical instructions and operating 
constants are then read in as data on punched cards. These are read in 
as variable names and are used in the various subroutines and are dis- 
cussed as part of the description of these subroutines. 
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Processing actually begins when MAIN reads in 4 additional in- 
structional constants on punched cards. These are given according to 
their variable name as 

NAME FUNCTION 

FILES - The number of the data file to be processed* 

ISTART - The record number within the file from which data 
will start being taken. 

ICHNO - The channel number of the magnetic tape to be 
processed. 

IRCEND - The record number at which processing is to stop* 

The values of these numerical constants enable the user to process any 
record or segment of records within a data file on any of the five 
channels. Each time a new file is to be processed these variables must 
be read in for the new file. The read statement for these variables is 
in a loop so that when processing of a file is completed the program re- 
turns to this statement to get instructions for processing the next file* 
After the last file has been processed the program is stopped by 
entering a negative number for the variable FILES. 

The program now moves the tape to the desired file and record by 
calling subroutine REDREC. When the first record to be processed is 
found, it is stored in a common array IBUF* MAIN then transfers the 
data in IBUF into a work array P. REDREC is called again to read the 
next record which is transferred by MAIN into an auxilary array AUX* 

This is done in order to have the next record available in an array 
since it is sometimes necessary for the program to "look" into the 
following record before processing in a record is completed. 

The record of data in array P is now processed by calling sub- 
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FORTRAN 

0001 

0002 


0003 

0004 

0005 

0006 

0007 

0008 

0009 

0010 
0011 
0012 

0013 

0014 

0015 

0016 

0017 

0018 

0019 

0020 
0021 
0022 
0023 


0024 

0025 

0026 

0027 

0028 

0029 

0030 

0031 

0032 

0033 

0034 


IV G LEVEL 1, MOD 4 (IAIN DATE = 70122 02/35/57 

INTEGER FILES FILCK 

COMMON P (1000) ,AUX( 1000) ,INDO(10) , NCNT (300) ,JOVF,JUNF, 

♦ BASE (2 f 10) , SIG (10) ,IBLOC (2,10) , KK , HH, XL90 , XLl 0 
*, DUMMY (2700) 

COMMON/ A A A/ 1 T BINT 

COMHON/AP/AilP (10000) ,NDATA 

COMMON/RLC/IRTN 

COMMON/RAT/IRCEND 

COM MON/GO/N 1 ,N3#N4,FM,IDO 

COMMON/NPTS8/N8 

COM EJQN/PR/I E (7,300) 

COM MON/O U U/N PNC H 
DIMENSION I Btl F (1000) 

DIMENSION IA (3) ,IB (5) 

CALL HID (1 1 , 1 2, IB) 

IP(I1) 155, 150,155 
155 PRINT 157, FILCK 

157 FORMAT ( 1 *♦* END OF FILE ♦♦♦•,18) 

150 IF (12) 158,161,158 

158 PRINT 159 

159 POR M AT ( 1 ♦♦♦ READ ERROR ♦♦♦•) 

161 CONTINUE 

NTAPE=IB (1) 

PRINT 160, IB 

160 FORMAT ( 1 TAPE NUMBER^ 1 , 1 A4 , 

1* 1DC AL- • ,1A4, 

1* BASE- 1 , 1A4 , 

It NUMBER OF CHANNELS' 1 , 1 A4 , 

1* NUMBER OF 5CANS=t,1A4) 

FILCK- 1 

READ (5, 9921) Ffl A X, DT, Ffi , N 1 , N3 , N4 , IDO , ISSB 
9921 FORMAT (3E10. 0,5110) ' 

READ (5, 3350) IPBINT,NPNCH 
3350 FORMAT (21 10) 

C READ INSTRUCTIONS AND TOLERANCES 

READ (5,36 0) IBTY P, I FLAG , NCI , NOSCAN, NOCRAN , LI , L2 

360 FOR fl AT (71 1 0) 

READ (5, 199 1) TOL 1 , TOL2, TOL3 

1991 FORMAT ( 3E 1 0 * 0) ^ ^ 

PRINT 361 ,IBTYP, IFLAG, NCI, NOSCAN , NOCK AN, TOL1 r TOL2,TOL3, LI r L2 

361 FORMAT ( ' TYPE BASE CALCULATION *************♦♦***********♦♦♦*♦♦♦ 

****** ***** IBTYP ***** *110/, 

*' TYPE STATISTICS CALLED FOR **********************♦***♦♦♦♦♦♦* 
****** IFLAG ***** «,I10/, 

*' NUMBER OF CLASS INTERVALS ***♦*************♦*♦*♦♦ *********** 
****** nci ***** * ,ri o/, 

*« NUMBER OF SCANS PER RECORD ****♦***********♦**♦***♦♦*♦♦♦♦♦+♦ 
****** HOSCAN ***** * ,110/, 

*» NUMBER OF CHANNELS ON TAPE ********************************* 
****** noCHAN ***** *,110/, 

*• TOLERANCE ON BASE LIMITS **+*****♦*********♦**************** 
****** TOL1 ***** • , Fl 0« 3/, 

** TOLERANCE ON SIGNAL LIMITS *********♦♦******♦♦**♦♦♦♦♦♦♦♦♦♦♦* 
****** TOL 2 ***** • , F 10.3/, 

*• NUMBER OF SIGNAL AND BASE POINTS REQUIRED PER CYCLE ******** 
****** XOL3 ***** * , FI 0 - 3/, 




FORTRAN IV G LEVEL 1 


HOD 4 


BAIN 


DATE = 70122 


02/35/57 


0035 

0036 


0037 

0038 

0039 

0040 

0041 

0042 

0043 

0044 

0045 

0046 

0047 

0048 

0049 

0050 

0051 

0052 


0053 

0054 

0055 

0056 

0057 
005B 

0059 

0060 
0061 

0062 

0063 

0064 

0065 


*» NUMBER OF PASSES REQUIRED TO ABORT CYCLE ******************* 

****** l 1 ***** '#110/, 

*• RUBBER OF EASE POINTS REQUIRED FOR SIGNAL POrNT SEARCH ***** 
****** 12 ***** * r I 10 ) 

PRINT 9 1 1 1,NPNCH,IPRINT,ISSfl r Ff1AX,DT,FH,Nl,N3 # N4,IDO 

9111 FOR MAT ( f PUNCHED OUTPUT FLAG *********************************** 
****** ***** NPNCH ***** *,110/# 

*• PRINT ERROR TABLE FLAG ****************** ****************** * 
****** IPRINT ***** »,I10/, 

** TYPE ANALYSIS ROUTINE WANTED ***♦ *************************** 
****** ISSB ***** *,110/, 

*» INFORMATION FOR FFT ROUTINE ************************ **+ **** * 
****** PM A Y ***** • ,F1 0. 3/,65X ***** DT ***** 

**,F10.3/ # 65X, '***♦• FM ***** ♦ r F 1 0 . 3/, 65X , • * **** Nl 

* ***** ' ,110/, 65X, '♦**** N 3 ***** ' ,110/, 65X, 1 **' 

**** N4 ***** ' ,110/, 65X, ' ***** IDO ***** Ml 

♦0/) 

C READ DATA FOR TAPE PROCESSING 

LFILE-0 

51 1 READ (5, 882) FILES, ISTART, ICHNO, IRCEND 
N 8 = 0 

N D AT A=0 . 0 

DO 8892 IJ=1,7 
DO 0892 J I x 1 , 300 
8892 IE (I J , JT) =0 

TRECNT=ISTART-1 
IEX IT=0 

IF (FILES-LFILE) 22,22,987 
987 TR EC=0 
22 CONTINUE 

LFILE=FILES 
882 FO R M AT { 4 1 3) 

PRINT 901, FILES, IC HNO , I 3T ART , IRCEND 
901 FORM AT ( 1 H 1 , 1 PROCESSING FILE ***********, 15/, 

*i CHANNEL ***************************** 1 , 15/ , 

*■ BEGIN AT RECORD ********************* 1 , 15/ , 

*• STOP PROCESSING AT RECORD *♦**♦**♦***', 15) 

I RTN-= 1 

IF (FILES) 991,9,9 

991 PRINT 992 

992 FORMAT (10X, f ****** PROGRAM STOP ******) 

STOP 

9 CALL KEDREC (FILES , I R EC , FI LC K , IBU F, TI M E , N , NOSC AN, NOCH A N, ICHNO) 

C MOVE TAPE TO DESIRED FILE 

IF (FILCK-FILES) 9, 887, 887 
C READ RECORD, STORE RECORD IN ARRAY IBUP 

887 CALL REDREC (FILES, I R EC , FTLCK , I BUF, TI ME , N , NOSCAN, NOCHAN, ICHNO) 

C MOVE TAPE TO DESIRED RECORD 

IF (I R EC- 1 START) € 87, 8 R 9, 8 89 
C STORE CONTENTS OF IBUF INTO MAIN ARRAY P 

889 DO 888 L= 1 , N 

888 P(L) =IBUF(L) 

C READ NEXT RECORD INTO IDUF 

CALL RED REC (FILES, IREC, FILCK, TBUF, TIME, N,N03 CAN, NOCHAN , ICHNO) 
c STORE CONTENTS OF TBUF INTO AUXILIARY ARRAY AD X 

DO S00 L= 1 , N 


assss* 



FORTRAN 17 G LEVEL 1, HOD 4 


MAIN 


DATE - 70122 


02/35/57 


0066 

0067 

0068 

0069 

0070 

0071 

0072 

0073 

0074 

0075 

0076 

0077 

0078 

0079 

0080 

0081 

0082 

0083 

0084 

0085 

0086 
0087 


500 AtJX (L) =1 8 OF <L) 

GO TO 501 

502 DO 503 L=1,R 

503 P (L) = AOX (L) 

CALL RED HEC<riLES,IREC,PILCK,IBUF, TINE, N,NQSCAN, NOCHAR, ICHNO) 
C STORE IBUF INTO AUX 

DO 504 L - 1 # N 

504 AUX (L) =IBUP(L) 

501 CALL EXT (N , L 1 , L2 , IREC NT , TOL1 , TOL2 , TOL3 ,NCI , I 3TTP) 

C STOP PROCESSING ON DESIRED RECORD 

IF ( IEXIT) 1000, 1000,507 
1000 IF (NDATA. EQ. 10000) GO TO 507 
661 IF (IR EC- IRC END) 502,221,221 
221 DO 2293 I=1,N 
2293 P(I)-AOXa) 

I EX IT= 1 
GO TO 501 

C PRINT ERROR TABLE 

507 CALL PRINT ( 1 ) 

IP (ISSfl- EQ. 2) GO TO 2514 
CALL FFT(N f DT,PJ1AX) 

IF (ISSB.EQ. 1) GO TO 511 
2514 CALL STAT.(NTAPE,ICHNO, FILES, NCI, IFLAG) 

GO TO 511 
END 
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routine EXT. When EXT completes the processing it returns program 
control to MAIN. The program checks the variable IEXIT to see if the 
record just processed is the last record in the file. If it is not, 
the program checks to see if the next record is the last. If this re- 
cord is not the last, the program calls in the next record and continues 
processing. However, if it Is the last record in the file to be pro- 
cessed, the program transfers the contents of the AUX array (which con- 
tains the last record) into work array P. IEXIT is set to 1 and EXT 
called to process the last record. When EXT returns program control to 
MAIN, IEXIT indicates that processing of this file has been completed. 

The program then calls subroutine PRINT to print the irregularity 
table. The analysis of the extracted signal continues by calling the 
statistical analysis subroutine STAT or the spectral analysis subroutine 
EFT. When the analysis is completed for this file of data, control is 
returned to MAIN which loops back to read the instructions for the next 
file to be processed. 

Subroutines 

A description of the subroutines called in the program is given in 
this section. The only subroutines not listed are FOURT and PLOT, since 
they wfere used in a package furnished by the computer department, A 
Fortran list is included after the description of each subroutine. 

SUBROUTINE RID (13, 14, IC) 

This subroutine reads the identification file which is the first 
file on the tape. This routine uses 3 subroutines that are especially 
written for the University of Alabama IBM- 360 Model 50 computer. The 
first, NTRAN, actually reads the tape and stores the data into an array. 
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The data is then transferred from this array and decoded for system 
compatibility using utility subroutines MOVE and TRNSL. This routine 
would probably require modification or rewriting if it were used on 
another machine* 

SUBROUTINE REDREC (FILES, IREC, FLICK, IBUF, TIME, N, NOSCAN, NOCHNO, 

ICHNO) 

Subroutine REDREC is called by MAIN to read and reformat the data 
from the magnetic tape* REDREC calls on the special utility subroutine, 
NTRAN to read the tape* Utility subroutines MOVE and TRNSL are called 
to convert the 7 track tape output into the byte system. REDREC must 
unpack from the multiplexed data the desired channel and convert the 
binary code to conventional base ten numbers. It must also keep up with 
the record and file number that it is reading. Provisions were made in 
REDREC to indicate read errors that might occur in NTRAN. This sub- 
routine would probably require modification if used on another machine* 


Argument Variables 

FILES - Number of files to be processed. 

IREC - Record number as counted by REDREC* 

FLICK - File number as counted by REDREC. 

IBUF - Array containing raw data* 

TIME - Not used. 

N - Number of data points per record (either 200 or 1000)* 

N0SCAN - Number, .of scans per record per channel (either 200 or 1000). 

NOCHAN - Number of channels multiplexed on the tape* 

Channel number to be processed. 


ICHNO 
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0001 


SUBROUTINE R TO (I 3 , 14 , IC) 

0002 


INTEGER 16 (5) ,IA (3) , BUP (501) ,FLCNT,TTB (3) , BLK, FILES , FILCK, 
*IB0F (1) ,IC (1) 

0003 


n=i 

0004 


n-o 

0005 


12-0 

0006 

5 

CALL NTRAN (1, 2,3 ,IA,K, 2,-2004 ,BUF,L,22) 

0007 

2 

. IF (K*1) 1,2,3 

0008 

3 

CALL MOVE {IB, 1,IA,1,4) 

0009 


CALL TRNSL (IB, 4, TTB) 

0010 


DATA TTB/*O123456709V 

0011 


data BLK /' v 

0012 


I B { 2) = BLK 

0013 


CALL MOVE (IB (2) , 1,IA,5,1) 

0014 


CALL TR N5L (IB (2) r 1 ,TTB) 

0015 


IB (3) = BLK 

0016 


CALL MOVE (IB (3) , 1,IA,6,1) 

0017 


CALL TRNSL (IB (3) r 1 , TTB) 

0018 


IB (4) “BLK 

0019 


CALL HOVE (IB{4) ,1,IA,7,2) 

0020 


CALL TRNSL (IB (4) ,2, TTB) 

0021 


CALL MOVE {IB (5) ,1,IA,9 ff 4) 

0022 


CALL TRNSX (IB (5) ,4, TTB) 

0023 


GO TO (6,7) , M 

0024 

6 

DO 14 1=1,5 

0025 

14 

IC(I)=IB(I) 

0026 


13 = 11 

0027 


14 = 12 

0026 


RETURN 

0029 

1 

IF (K.EQ.-2) GO TO 4 

0030 


12=1 

0031 


CALL NTH AN (1,22) 

0032 


GO TO (6,7) 

0033 

4 

11 = 1 

0034 


CALL NTRAN (1,22) 

0035 


GO TO (6,7) ,K 

0036 


ENTRY REDR EC (FILES, I R EC, PILCK , I 8UF,TI HE, N, NOSCAN, NOCHAN, ICHNO) 

0037 


tt = 2 

0038 

13 

IREC=IR£C+ 1 

0039 


CALL NTRAN (1,22) 

0040 

9 

IF (L* 1)8,9,10 

0041 

10 

TIM E=BUF ( 1 ) 

0042 


N= (L— 4) / (2*NOCHAN) 

0043 


K=3+2*ICHNO 

0044 


DO 11 1=1, N 

0045 


I BUF (I) =0 

0046 


J = 0 

0047 


CALL MOVE (3, 4, BUF (1) ,K,1> 

0048 


CALL MOV E (IBU F (I) ,4,BUF(1> ,K*1 f 1) 

0049 


IBUF(I)*IBUF (I) + 64* J 

0050 


IF (IBOF (I) .GE. 1024) I BOF (I) = 1 02 4- IBUF (I) 

005** 

11 

K = K + 2*N0CH AN 

0052 


CALL NTRAN (1, 2, -2004, BUF, L) 

0053 


RETURN 

0054 

8 

IF (L.EQ.-2) GO TO 12 

0055 


PRINT 100 


ORIGINAL PAGE IS 
OF POOR QUALITY 
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0056 

100 

PORflAT (»0****** READ ERROR 

0057 


CALL NTRAN (1, 22) 

0058 


CALL NTRAN ( 1 , 2,-2004 , BU F ,L) 

0059 


GO TO 13 

0060 

12 

CONTINUE 

0061 


FILCK“FILCK* 1 

0062 


CALL NTR AN (1,22) 

0063 


GO TO 5 

0064 

7 

CONTINUE 

0065 


IREC-0 

0066 


GO TO 13 

0067 


END 


ORIGINAL PAGE tS 
OF POOR QUALM 
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SUBROUTINE EXT (NK, Ll, L2 , IREC, TOL1, TOL2, TOL3, NCI, IBTYP) 

This subroutine is by far the most complex routine in the program* 
It extracts the intensity amplitudes from the chopped data. A written 
explanation of EXT will not be given due to its complexity. Included 
Instead is a flow diagram. It is hoped that the interested reader can 
use the description given in the text along with the flow chart to 
understand the operation of this routine. As an additional aid, the 
important variable names are given a brief description. 


Argument Variables 


NK - Number of data points per record. 

LI - Number of sequential searches for base and signal points 

allowed before cycle is aborted. 

L2 - Number of base points required for projected signal point 

( search. 

IREC - Record number being processed. 

T0L1 - Sets tolerance on base limits for base point selection. 

T0L2 - Sets tolerance on signal limits for signal point selection. 

T0L3 - Number of signal and base points required for a normal cycle. 


NCI - 


Number of class intervals. See subroutine HIST. 


IBTYP - Type base calculation. See subroutine HIST. 


Common Block Variables 

P - Work array containing record of data points being processed. 

AUX - Auxiliary array containing the next record to be processed. 

INDO - Array containing position in the data array from which the 
signal points came. 

JOVF - Number of overflows in HIST. Variable is initialized in EXT. 

JUNF - Number of underflows in HIST. Variable is initialized in EXT. 

BASE - Array containing both groups of base points. 
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SIG - Array containing signal points • 

I BLOC - Array containing position in the data array P from which 
each group of base points came. 

L, M - When EXT calls subroutine LIMIT, it gives it the initial 
position L and the final position M LIMIT is to scan in 
the data array. 

XL90 - This is the result of calling LIMIT and is the criteria 
for selecting base points. 

XL10 - Also the result of LIMIT and is the criteria for selecting 
signal points. 

Labeled Common Variables 

IPRINT - If IPRINT is other than 0, EXT will print the record if any 
irregularities occur while the record is being processed. 

If IPRINT = 0 it will avoid printing, 

IRTN - Place keeper for subroutine EXT. The subroutine may be in 
any part of its cycle when it completes a record since data 
is continuous from record to record. IRTN is set to a number 
corresponding to the exit point in the routine when it re- 
turns to the MAIN routine for a new record. When subroutine 
EXT is recalled, a computed GO TO statement keyed to IRTN 
returns control to the phase EXT was previously in. 

IE - Array passed to subroutine print which contains the errors 

accumulated for each record. 

IBCNT - Array containing the number of base points in each group. 

ISCT - Number. of signal points. 


Important Internal Variables 


IRUN - IRUN = 1 indicates routine in start up cycle. Converse for 
IRUN = 0. 


INDX - The value of this variable indicates the position (1-200) in 
the record being processed. 

IDINX - The number of positions subroutine FILL must place zeros due 
to irregularities which cause data points to be skipped. 

NP - Count of unsuccessful passes through signal and base search 

cycle. 


IBC - 


Indicates which group of base points are being searched for. 
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LC - Indicates loop condition* LC - 0 means routine is the 

base point search phase; LC == 1 indicates signal point 
search* 

INP - Counts the times the error recycle phase is entered* 

SUBROUTINES HIST (BA, IDUM, NCI, IBTYP , I RUN) 

This routine is called by EXT to compute the amplitude of the 
chopped wave and constructs a histogram with the results* The routine 
is designed to use numbers between 0 - 1000 but can be easily modified 
to handle a larger range. The histogram is stored in a common array 
to be used by the statistical analysis subroutine STAT. 


Argument Variables 


BA - Average of the base points . 

NCI - Number of class intervals for histogram. 

IBTYP - Determines the method to be used to calculate the amplitude. 
If IBTYP = 0, amplitudes are calculated by taking the 
difference between the signal points and the average of the 
base points in both groups. If IBTYP » 1, calculation will 
be the difference between the first signal point and the last 
base point of group 1, and the difference between the last 
signal point and the first base point of group 2. If 
IBTYP = 3, calculation is performed only on the last base 
point of group 1 and the first signal point. 


IRUN - If equal to 1, indicates that EXT ; is in.' its "first -cycle. 


Common Block Variables 

NCNT - Array containing frequency for each class interval. 

JOVF - Number of overflows or excessively large numbers resulting 
from classifications of amplitudes. 

JUNF - Number of underflows or very small numbers resulting from 
classification of amplitudes. 

BASE - Array containing two groups of base points. 

SIG - Array containing signal points. 
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FORTRAN 

0001 

0,002 

0003 

0004 

0005 

0006 
0007 
OOOB 

0009 

0010 
0011 
0012 

0013 

0014 

0015 

0016 

0017 

0018 

0019 

0020 
0021 
0022 

0023 

0024 

0025 

0026 

0027 

0028 

0029 

0030 

0031 

0032 

0033 

0034 

0035 

0036 

0037 

0038 

0039 

0040 

0041 

0042 

0043 

0044 

0045 

0046 

0047 

0048 

0049 

0050 

0051 

0052 


IV G LEVEL 1, HOD 4 EXT DATE = 70122 

SUBROUTINE EXT (NK,Ll , L2, IREC,TOL1 # TOL2 ,TOL3 , NCI, IBTYP) 
COHflON P (1000) ,AUX( 1000) ,INDO(10) ,NCNT(300) ,JOVP,JUNF, 
♦ BASE (2, 10) ,SIG (10) ,IBLOC (2,10) ,L, M, XL90,XL 10 
COHHON/AAA/IPRINT 
COMflQN/BLC/IETN 
COM HON/PR/IE (7,300) 

COM NON/CBA/I BCNT (2) ,ISCT 

IREC=IRECM 

IST=IST ♦ 1 

INP-O 

INO=1 

I NDX= 1 

GOTO (150, 11,12,150,776) # IFTM 
150 IRU N= 1 
I RP = 0 
1ST = 0 
L - 1 
H = 20 

IF (IRTN.EQ.4) GO TO 209 
DO 1050 1=1,300 
1050 NCNT(I)=0 

JOV?=0 
JUNF=0 
NDATA=0 

209 CALL LIBIT (TGL 1 , TOL2 , NK) 

C SEARCH FOR FIRST BASE POINT- HOLD VALUE OF INDX 

776 DO 1 1= L, M 

I F ( I. GT . N K) GO TO 2000 
IF (P (I) -XL90) 1,3,3 
3 INDX=I 

GO TO 4 

I CONTINUE 

IE (1,IREC)=IE (1 ,IREC) *1 
IRTN=4 

C RETURN TO DRIVER FOR NEXT RECORD 

RETURN 
2000 L=1 

N=M-NK 
IRTN-5 
RETURN 
4 H=INDX+15 
L=INDX 

CALL LIWIT (TOL 1 , TOL2 , NK) 

I DI NX=I ADS (INDX-INO) 

CALL FILL (I DINX) 

INO=INDX 
IBC= 1 

I BCNT (1) -0 
I BCNT (2 ) =0 
ISCT=0 
10 LC = 0 
NP= 0 

C BASE POINT SEARCH 

II IF (P (INDX) -XL90) 12,13,13 

13 IBCNT(IDC)=IBCNT (IBC) +1 

I = T BCNT (IBC) 
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FORTRAN 

IV G LEVEL 

1, MOD 4 EXT DATE 


C 

STORE BASE POINT 

0053 


BASE (IBC, I) = F (I NDX) 

0054 


I BLOC (IBC, I) =INDX 

0055 


IF<IBCNT{IBC) -10) 14,401,401 

0056 

401 

IRTN*4 

0057 


I E (2,IR EC) -IE (2, IREC) +1 

0058 


IDI NX- NK- I NDX 

0059 


CALL PILL (IDINX) 

0060 


RETURN 

0061 

14 

INDX*INDX ♦ 1 

0062 


INO-INDX 

0063 


IF {INDX- NK) 11, 11,158 

0064 

158 

IRTN-2 


C 

NORMAL RETURN TO DRIVER FOR NEXT RECORD 

0065 


RETORN 


c 

SIGNAL POINT SEARCH 

0066 

12 

IF {P (IN OX) -XL 10) 20,20,15 

0067 

15 

I F ( LC) 16, 16,21 

0068 

16 

NP=NP*1 


C 

CHECK NUMBER OF UNSUCESSFUL PASSES 

0069 


IF(NP-LI) 1944, 1944, 200 

0070 

1944 

IF (ISCT.EQ.O) GO TO 14 

0071 


GO TO 70 \ 

0072 

21 

IF ( P (INDX) -XL90) 22,3 0,3 0 

0073 

30 

IF (IRON) 10,10,31 


C 

SET BASE COUNT INDX CONDITION 

0074 

31 

IF(IBC-I) 34,33,34 

0075 

33 

IBC -2 

0076 


GO TO 10 

0077 

34 

I BC- 1 

0078 

20 

IF (IRON) 23,23,24 


C 

CHECK LOOP CONDITION 

0079 

2 3 

IF (LC) 40,40,29 

0080 

24 

IF(IBC-I) 2 5,25,40 

0081 

2S 

LC= 1 

0082 


NP" 0 

0083 

29 

ISCT=ISCT ♦ 1 


C 

STORE SIGNAL POINT 

0084 


5IG (ISCT) -P (INDX) 

0005 


INDO (ISCT) = INDX 

0086 


IF (ISCT* 10) 26,501,501 

0087 

501 

IBTN-4 

0088 


IE (3, IREC) -IE (3, IREC) *1 

0089 


idinx^nk-indx 

0090 


CALI FILL (IDINX) 

0091 


RETURN 

0092 

26 

INDX-INDX+1 

0093 


INO-INDX 

0094 


IF(INDX-NK) 12,12,162 

0095 

162 

IRTN=3 


C 

NORMAL RETURN TO DRIVER FOR NEXT RECORD 

0096 


RETURN 

0097 

22 

NP=NP+1 

0098 


IF* (NP-L1) 26,26,208 


C 

CHECK NUMBER OF BASE AND SIGNAL POINTS 

0099 

40 

IF(IBCNT(1) pLT.TOL3) GO TO 222 


«SS*SS“ 
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FORTRAN IV G LEVEL 1, HOD 4 


EXT 


DATE = 70122 


0100 

0101 

0102 

0103 

0104 

0105 

0106 

0107 

0108 

0109 

0110 
0111 
0112 

0113 

0114 

0115 

0116 

0117 

0118 

0119 

0120 
0121 
0122 

0123 

0124 

0125 

0126 

0127 

0128 

0129 

0130 

0131 

0132 

0133 

0134 

0135 

0136 

0137 
013B 

0139 

0140 

0141 

0142 

0143 

0144 

0145 

0146 

0147 

0148 

0149 

0150 


IP {IBCNT (2) , LT.TOL3) GO TO 222 
IF (ISCT.LT. TOL3) GO TO 224 
GO TO 100 

222 IE(4,IREC) -IE (4, IREC) *1 
GO TO 555 

224 IE {5, IREC) =IE (5, IREC) ♦! 

GO TO 555 

C ENTER ERROR RECYCLE 

208 IE (6,TREC)=IE (6, XREC) +1 
555 L*TNDX 

pi = INDX*20 


IRU N = 1 

IF (IPRINT. EQ.O) GO TO 742 
PRINT 666, IREC, INDX 

666 FORMAT (IX, • IflEC=* ,110, 10X, * INBX- * #110) 
IF (IREC.EQ. IRP) GO TO 742 
PRINT 333, <P(I) 

333 FORMAT ( IX ,/ r (IX, 1 OP 10. 3) ) 

IRP - IR EC 
742 CONTINUE 
I NP-INP * 1 

IF (INP-5) 209,921,921 
921 IE (7, IRECJ "IE (7, IREC) 


I RTN - 4 
IDINX=NK-INOX 
CALI FILL (IDI NX) 

RETURN 

POINT SEARCH CYCLE COMPLETE 


100 SUM 1 = 0 . 0 
SUM 2 = 0,0 


C 


101 

102 


43 

45 


1=1 BCNT (1) 

PREPROSSING FOR HISTOGRAM FOLLOWS 

DO 101 J=1,I 

SUM 1 = S0M 1+ BASE (1,3) 

I=IBCNT (2) 

DO 102 J=1,I 

SUM 2-SUM 2 *BA5E(2,J) 

BA= (SUM US UM2)/ (IBCNT (1) ♦IBCNT (2)) 
CALL AMPX (IRUN) 

CALL HIST (BA,NK,NCI,IBTYP,IRUN> 

IF ( IBC- 1) 44,45,44 

IBC-2 

GO TO 46 


44 I BC= 1 

46 I SCT“0 

IBCNT (IBC) “0 
LC- 1 


NP=0 


C 

70 


IRUN=0 

IN D X 2~ INDX ♦ 9 

L-INDX 

M=I NDX2 

CALL LIMIT (TOL1,TOL2, NK) 

GO TO 12 

LOOK INTO NEXT RECORD FOR SIGNAL POINT 
IF (IBCNT (IBC) -L2) 14,72,72 
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FOBTRAR IV G LEVEL 1, BOD 4 EXT DATE = 70122 

0152 72 L=I BOX 

0153 H=ISDX ♦ 9 

0154 CALL LIBIT (TOLl ,TOL2, HK) 

0155 IF(P(INDX|-XL10) 20,20,14 

0156 BSD 
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FORTRAN IV G LEVEL 1 , MOD 4 HIST 


DATE = 70122 


0001 

0002 

0003 

0004 

0005 

0006 

0007 

0008 

0009 

0010 
0011 
0012 

0013 

0014 

0015 

0016 

0017 

0018 

0019 

0020 
0021 
0022 

0023 

0024 

0025 

0026 

0027 

0028 

0029 

0030 

0031 

0032 

0033 

0034 

0035 


SUBROUTINE HIST (BA,TDUH , NCI , IBTYP, IRON) 

COMMON P (1000) ,AUX (1000) , INDO (10) , NCNT (30 0) ,JOVF, JUNF, 
♦ BASE (2, 10) ,SIG (10) , I8LOC (2 , 1 0) , KK , MW , XL9 0 , XL 1 0 
DIMENSION 1(40), X (40) , LOC (40) , A(15) f V(15), B (200) 
COM ttON/CBA/IECNT (2) tf N 
C COM BINE BASE 1 AND BASE 2 

IF ( IBTYP) 100, 100, 200 
200 IF(IRUN.EQ. 1) MMM = 1 
KKK - 1 
LLL= 1 

IF ( WMM- EQ . 1) LLL-2 
IF ( MMM . EQ . 0) KKK= 2 , 

100 IJ=N 

DO 1 I- 1 f N 

IF (IBTYP* EQ* 0} GO TO 202 

400 IF (I- 1) 402,401,402 

401 HTEMP=IBCNT (KKK) 

BA-BASE (KKK,NTEKP) 

GO TO 202 

402 IF(I-IJ) 1,403,1 

403 BA^BASE (LLL, 1) 

202 J = B A-SIG (I) 

J= (NCI*J)/1000 *1 

IF(J-I) 2,3,3 

2 JUNF=JUNF+1 
GO TO 1 

3 IF (J-300) 5,5,4 

4 JOV F« JOV P+1 
GO TO 1 

5 NDATA=NDATA+ 1 
NCNT (J) - NCNT (J) ♦ 1 

IF (IBTYP. EQ-3) GO TO 50 
1 CONTINUE 

IF (flrtPl.EQ.O) MtlM-1 
IF (MMN.EQ, 1) MMH^O 
50 RETURN 
END 
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SUBROUTINE STAT (NT APE , NCH, NFILE, NCI, I FLAG) 

STAT performs the statistical analysis by using the intensity 
histogram constructed by HIST. A log-amplitude histogram is generated 
from the intensity histogram to perform log-normal tests. The mean, 
standard deviation, skewness and kurtosis are calculated for both the 
intensity and log-amplitude data. In addition, the cumulative probability 
is calculated and a chi-square test made on both sets of data. 



Common Block Variables 

NCNT - Array contains histogram. 

JOF - Number of overflows. 

JUF - Number of underflows. 


SUBROUTINE LIMIT (T0L1, T0L2, NK) 

This routine is called by subroutine EXT to calculate the criteria 
for determining if a data point is a base point or a signal point or 
neither. LIMIT has the capability of looking into the next record if it 
is called near the end of the record being processed. 
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0021 

0022 

0023 

0024 

0025 

0026 

0027 

0028 

0029 
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IV G LEVEL 1, MOD 4 


STAT 


DATE = 70122 02/35/57 


C 

C 

c 

c 


c 


c 


c 


c 


c 


c 


c 


SUBROUTINE STAT (NT APE, NCH, N PILE , NCI , X FLAG) 

IFLAG = I NO CHI SQUARE TEST 

2 CHI SQUARE TEST ON NORMAL DISTRIBUTION 
.3 CHI SQUARE TEST ON LOG NORiILL DISTRIBUTION 
4 CHI SQUARE TEST ON BOTH 

COM HON P(1000),AUX(1000),INDC(10),NCNT1(300) , JOF, JUF, BASE (2 # 10) , 
♦PIG (10) , I B LOC (2,10) ,KK,HM, XL90, XL 10 
COB HON/UUU/N PNCH 

DIMENSION I (300) ,XLN (300) ,Q (4) ,R (4) 

DIMENSION NCNT (3 00) 

DO 60 1*1,300 
60 NCNT (I) =NCNT 1 (I) 

C- 1 000/NCI 

FIND THE HIGHEST and LOWEST CLASS INTERVAL 
DO 1 1=1, NCI 
IF (NCNT (I) ) 1,1,2 

2 ILO= I 
GO TO 3 

1 CONTINUE 

3 DO 4 1= ILO, NCI 
IF (NCNT (I) > 4,4,5 

5 I HI = I 

4 CONTINUE, 

FIND NUHBER OF POINTS AND FLOAT NCNT 
N = 0 

DO 6 I=ILO, IHI 
Y (I) =NCNT (I) 

6 N=N+NCNT ( I) 

PRINT HEADINGS 

WRITE (6, 101) NT APE, N CH , NFILE 
WRITE (6,102) 

DEFALT DUE TO TOO FEW CLASS INTERVALS 
N N= IHI-ILO 
IF(NN-IO) 50,50,51 

50 WRITE (6,110) 

RETURN 

FIND AVERAGE AMPLITUDE 

51 X N = N 
AVE=0.00 

DO 8 I=ILO,IHT 

XI*I 

8 AV E = AVE + Y (I) * ( XI~ 0. 5) *C 

A V E= A VE/X N _ . 

COMPOTE CUMULATIVE PROS ABILITIEC AND LOG AMPLITUDES 

3UH=0. 0 0 

DO 10 I=ILO, IHI 

XI=I 

XT = (XI-0. 5) ♦C 

XLN (I) = 0- 5* ALOG (XI/AVE) 

SUM=SUM*Y (I) 

CP=SUM/XN 

10 WRITE (6, 103) XI, XLN (I) , Y (I) ,CP 
WRITE(6,111) JOF, JUF 
COMPUTE MOMENTS ABOUT THE MEAN 
XL A= 0- 0 0 
DO 20 I = I LO, IHI 
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STAT 


DATE = 70122 


02/35/57 


0045 

0046 

0047 

0048 

0049 

0050 

0051 

0052 

0053 

0054 

0055 

0056 

0057 

0058 
0054 
0060 
0061 
0062 

0063 

0064 

0065 

0066 
0067 
006B 

0069 

0070 

0071 

0072 


0073 

0074 

0075 

0076 

0077 

0078 

0079 

0080 
0081 
0002 

0083 

0084 

0085 

0086 
0087 


0088 

0089 


20 XLA=XLA*Y (I) *XLN (I) 

XL A=XLA/XN 

DO 21 1-2,4 
0 | I ) =0.00 

21 R (I) =0.00 

DO 22 I=ILO,IHI 
XI=I 

DO 22 J = 2,4 

Q (J> =Q (0) *Y (I) * ( {XI-0. 5) +C- AVE) **J 

22 H ( J) =R ( J ) * Y (I)* (X LN <I)-XLA)**J 
DO 23 J = 2, 4 

Q(J)=Q(J)/XN 

23 R <J) =R(J)/XN 
NT=IHI-IL0+1 
WRITE (6, 10 4) NT 
SIG^SQRT (Q (2) ) 

SIGt=SQRT(R (2)) 

SK EW = 0. 5*Q (3) / (SIG**3) 

SKL = 0. 5*R (3)/(5IGL**3) 

XKU R- ( (0(4) / (Q(2) *+2) )-3.0)/2.0 
XKURL = ( <R (4)/ (ft (2) **2) >~3.0)/2.0 
C PRINT MOMENTS 

WRITE (6,105) AVE,SIG, SKEW , XKUR ,XLA,SIGL,SKL, XKURL, N 
I F ( NPNCH .EQ. 0) GO TO 810 

PUNCH 800, NT APE, WCH , NFILE, A VE, SIG,SIGL, XLA 
800 FORMAT (A 4,12,14, 4E 14. 4) 

810 CONTINUE 

GO TO (31,32,42,32) , IFLAG 

31 RETURN 

C NO CHI SQUARE TEST REQUESTED 

C 

C CHI SQUARE TEST 

32 CALL CHI (CSQ, Y , ILO , I HI, C , NU5E , A VE , XL A , SIG ,XN,0,SIG) 

X N = CSQ 

WRI T E (6 , 1 06) CSQ,NUSE 
C PRINT CHI SQUARE NORHAL 

GO TO (31,31,42,42) , IFLAG 

42 CALL CHI (CSQ, Y , ILO , IHI ,C , NUSE , A VE , XL A , SIG ,XN,1,SIGL) 

AXL N=CSQ 

WRITE{6, 106) CSQ , N USE 
C PRINT CHI SQUARE LOG NORMAL 

PUNCH 999, NTAPE, NCH , NFI LE , SIGL, SKEW, SKL,XN,AXLN 
999 FORMAT (A4, 1 1,12, 5E 13. 8) 

RETURN 

101 FORMAT (• 1 » ,5X, 'TAPE NUMBER *, A4, 5X, 'TRACK*, 13, 5X,' FILE* ,13) 

102 FORMAT ( * 0 1 ,16X, * AMPLITUDE 1 ,10X, 1 LOG AM P L ITUDE ' , 1 2X, » COU NT » , 

1 8X, * CUMULATIVE PROBABILITY'/) 

103 FORMAT (7X,4E21. 6) 

104 FORMAT ( ' 0 * , 10X, 'NUMBER OF CLASS INTERVALS =*,I6) 

105 FORMAT ('0* ,17X, 'AVERAGE' , 9X , ' ST A ND ARD DEV I ATI ON • , 8X , ' SK EW NESS * , 
*13X, 'KURTOSIS */7X, 4 E2 1 . 6/7X , 4E2 1 .6/'0* ,10X, • NUMBER OF DATA POINTS' 
*,H0) 

106 FORMAT ('O' , 10X, 'CHI SQU A RE= 1 , El 4 . 6/ 11X ' NUMBER OF CLASS INTERVALS 
1 USED = ' , 15) 

110 F0RMAT('0' ,5X,' TOO FEW CLASS INTERVALS '/ 

1 *0*,5X,' EXECUTION OF STATISTICS CALCULATION SUSPENDED') 
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FORTRAN IV G LEVEL 1 , NOD 4 STAT DATE = 70122 

0090 111 FORNATrO* ,5X ff 'NUNBER OF OVERFLOWS 1 ,16/ 

*6X,'NUNBER OF UNDERFLOWS • f IS) 

0091 EKD 




02/35/57 
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Common Block Variables 

KK - Gives the position (value of INDX) in the record at the 

time LIMIT is called. 

MM - This variable is the sum of KK and the number of points 

LIMIT is to scan. 

XL90 - The resulting criteria for base point selection. 

XL10 - The resulting criteria for signal point selection. 

SUBROUTINE AMPX (IRUN) 

This routine takes the signal and base points extracted by EXT and 
computes the amplitude of the square wave for the spectral analysis. The 
amplitudes are calculated by taking the difference between the last base 
point in the first group and the first signal point, and the difference 
between the last signal point and the first base point in the second group. 
This produces two signal points per group. The points are stored in 
array AMP for use by the spectral analysis routines. 

Argument Variables 

IRUN - Indicates if EXT in startup cycle. 

Common Block Variables 

BASE - Array containing both groups of base points. 

SIG - Array containing signal points. 

Labeled Common Variables 

IBCNT - Array containing number of base points in each group. 

N - Number of signal points 

SUBROUTINE CHI (CSQ, Yl, ILO, IHI, C, NUSE, AVE, XL A, SD, XN, NTYP , SX) 

This routine is called by the statistics subroutine to perform the 
chi-square test for normal and log-normal distributions. 
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0001 


SOD ROUT IN B LIMIT (TOL1 ,TOL2,NK) 

0002 


COMMON P<1000) ,AUX (1000) , IN DO (10) ,NCNT (300) r JOVF, JTJNF, 
♦ BASE (2,10) , SIG ( 1 0) , I BLOC (2, 10) , K K , MM , XL90 , X L 1 0 

0003 


IK=KK 

0004 


IM = MM 

0005 


ID^Q 

0006 


IF (IK-NK) 502,501,666 

0007 

666 

PRINT 3000 

0008 

3000 

FORMAT ( * IK GREATER THAN NK ') 

0009 

501 

AMAX^P (IK) 

0010 


AM I N- AM AX 

0011 


ID"= IM-NK 

0012 


GO TO 381 

0013 

502 

IF(IM-NK) 360,360,361 

0014 

361 

ID=IM-NK 

0015 


IH-NK 

0016 

360 

AM AX-P (IK) 

0017 


AMIN^AMAX 

0018 


DO 350 J=IK,IM 

0019 


I F ( AM AX— P ( J) ) 301,302,302 

0020 

301 

A M AX = P (J) 

002 1 

302 

TP (AM IN- P (J) ) 350,303,303 

0022 

303 

AMI N=P (J), 

0023 

350 

CONTINUE 

0024 


TF(ID) 380,300,381 

0025 

301 

AM A XP = AUX (1) 

0026 


AMINP^AMAXP 

0027 


DO 650 J= 1 f ID 

0028 


IF (AMAXP-AUX (J) ) 601,602,602 

0029 

601 

AMAXP-AUX (J) 

0030 

602 

IF (AMINP-AUX (J) ) 650,60 3,603 

0031 

603 

AM I NP- AUX (J) 

0032 

650 

CONTINUE 

0033 


I F ( AM AX- AM AXP) 700,701,701 

0034 

700 

A M A X- AM AX P 

0035 

701 

IF (AMIN-AMINP) 380,380,703 

0036 

703 

AMIN=AMINP 

0037 

380 

A -A H AX' AMIN 

0030 


XL90-AMAX-TGL1*A 

0039 


X L 1 0 = AMI N + TOL 2* A 

0040 


RETURN 

0041 


END 
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0001 


SUBROUTINE AMPX (IRUNJ 

0002 


COM WON P( 1000) , A U X (1000) ,INDQ(10) ,NCNT{300) ,JOVP,JUNF, 
♦ BASE (2, 10) ,SIG(10) , IBLOC (2,10) , KK f MM , XL90 , XL1 0 

0003 


COMMON/CBA/IBCNT (2) ,N 

0004 


COMMON/AP/AMP (10000) , NBA T A 

0005 


IP (IFUN. EQ. 1) MrtM=1 

0006 


KKK = 1 

0007 


LLL- 1 

0008 


IP (MMM.EQ. 1) LLL-2 

0009 


IF (MHM. EQ. 0) KKK = 2 

0010 


I J = N 

0011 


DO 1 1=1, N 

0012 

400 

IF (1-1) 402,401,402 

0013 

401 

NT EMP=I BCNT (EKE) 

0014 


B A= 8 AS E (KKK , NTEM P) 

0015 


GO TO 502 

0016 

402 

IF(I-IJ) 1,403,1 

0017 

403 

BA-BASE (L1L, 1) 

0018 

502 

NDATA~NDATA* 1 

0019 


AMP (NDATA) =BA-SIG (I) 

0020 

1 

CONTINUE 

0021 


IP(MMM.EQ.G) MMM= 1 

0022 


I F ( M MM . EQ . 1) MMM^O 

0023 


RETURN 

0024 


END 
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CHI 


DATE * 70122 


02/35/57 


0001 

0002 

0003 

0004 

0005 

0006 

0007 

0008 

0009 

0010 
0011 
0012 

0013 

0014 

0015 

0016 

0017 

0018 

0019 

0020 
0021 
0022 

0023 

0024 

0025 

0026 

0027 

0028 

0029 

0030 

0031 

0032 

0033 
0014 
00.15 

0036 

0037 
0030 
0039 


SUBROUTINE CHI (CSQ, Y1,ILO,IHI,C, NUSE r A VS , XL A , S D r XN, NTYP, SX) 
DIMENSION 7(300, 3) . Y 1 (300) 

DD 1 1=1,300 

1 Y (1 , 1) -0. 00 
NUS£=0 

KM = (ILO+IHI)/2 
J=1 

Y (3,2) = A V E~ 1 0 . 0 *SD 

C GROOP CLASS INTERVALS ON LOW END 

DO 2 1= tLO.KR 

Y (3, 1)=Y1 (I) + Y (3,1) 

IF (Y (J, 1) -5.0) 2,2, 3 

3 Y(3,3)=C*t 

NUSE - NUSE ♦ 1 
J = J ♦ 1 

Y (3,2) = C*I 

2 CONTINUE 

C GROOP CLASS INTERVALS ON HIGH SIDE 

I = IHI 

' Y (3,3) = A VE + 10.0+SD 
6 IP(I-KM) 10,10,4 

4 Y (3, 1) = Y (J,1) + Y1 (I) 

I F ( Y (3 , 1).-5. 0) 11,11,5 

5 Y(J,2) = 

NUSE = NUSE ♦ 1 
3 = 3+1 

Y (J , 3) = C*(I-1) 

111=1-1 

GO TO 6 

C COMPUTE THEORITICAL PROBABILITY 

10 CSQ = 0.00 

DO 30 1=1, NUSE 
XLL = T (1,2) 

XUL = Y (1,3) 

24 CALL SIMP ( FTH , XLL , XU L , 2 1 , NT7P , A VE , SX , XL A) 

C COMPUTE CHI SQUARE 

IF (FTH) 31,31,30 
31 WRITE (6, 100) I 

100 FOR M AT ( ■ 0 ZERO VALUE OF THEORITICAL PROBABILITY IN 1 , 15,* TH 
1 INTERVAL’/ 6X, 1 EXECUTION OF CHI SQUARE TEST DISCONTINUED 1 ) 
RETURN 

30 CSQ=CSQ + ( (Y (I, 1)-XN*FTH) ** 2 ) / (X N* PTH) 

RETURN 

END 


ORIGINAL PAGE IS 
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Argument Variables 

CSQ - Result of chi-square test. 

Y1 - Array containing histogram. 

ILO - Lowest class interval in histogram. 

IHI - Highest class interval in histogram. 

C - Width of class mark. 

NUSE - Number of class intervals used by the chi-square routine. 

AVE - Mean value of amplitudes. 

XLA - Mean value of log-amplitudes . 

SD - Standard deviation of amplitudes. 

XN - Number of data points. 

NTYP - Determines if chi-square test will be run for normal or log- 
normal test or both. 

SX - Log standard deviation. 

SUBROUTINE FFT (DT, FMAX) 

This routine is called by the main program to coordinate the per- 
formance of the spectral analysis. Subroutines FOURT and PLOT are called 
to perform the Fourier transform and to plot the results. FFT will have 
PLOT plot directly from the calculated spectral data array or it will have 
it plot the average of a designated number of points in the array. This 
feature was incorporated to smooth out random variations. 

Argument Variables 

DT - Time between data samples. 

FMAX - Maximum frequency to be used. 

Labeled Common Variables 

AMP - Array containing the time domain signal. This array is 

passed from AMPX. 
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0001 


SUBROUTINE FFT (N , DT r FH AX) 

0002 


DIMENSION WORK (2 500) 

0003 


COM WON/ A P/AMP {10000) , N DAT A 

0004 


COM MON/GO/ Nl,N3,N4,FJ1 r IDO 

0005 


C0MJ10N/NPTS8/N8 

0006 


DIMENSION NN<2) 

0007 


PRINT 555, NDAT A 

0008 


. IF (NDATA-8192) 25, 56 , 56 

0009 

25 

NDIFF^8192-NDATA 

0010 


CALL FILL { NDI FF) 

0011 

56 

PRINT 81 r N0 

0012 

81 

FORKAT{10X, f NUMBER OF 2 EROS ADDED = *,T20//) 

0013 


N DAT A=8 19 2 

0014 

555 

FORMAT ( 10 X , * NUMBER OF DATA POINTS EXTRACTED -*,120,///) 

0015 


N~ N DAT ft 

0016 

99 

FORMAT (/, * IFREQ(HZ) MAGNITUDE* ) 

0017 


NN (1)=N 

0018 


DF = 1 ■ 0/ (N*DT) 

0019 


N 2 - (PMAX/DF) *2 

0020 


IP (N2. GT- N) N 2 = N 

0021 


CALL FOURT (AMP, NN, 1,-1, 0, WORK, 2500) 

0022 


DO 5 J-2,N,2 

0023 


X = A HP (J- 1) 

0024 


T” A MP (J) 

0025 

5 

AMP (J-1) -SQRT (X*X+Y*Y}/N 

0026 


WRITE (6,99) 

0027 


IF ( IDO. EQ* 1) GO TO 50 

0028 


CALL PLOT (AMP, DF,PF,N1 ,N2,N3) 

0029 

50 

F = D F 

0030 


IF (IDO.EQ.3) GO TO 51 

0031 


PRINT 100 

0032 

100 

FORMAT (/,* FREQ (HZ) MAGNITUDE 1 ) 

0033 


SUM -0 « 0 

0034 


J-1 

0035 


NC= 0 

0036 


DO 10 I=N1,N2,N3 

0037 


NC = NCM 

0038 


SU K” AM P (I) + SUH 

0039 


IF(NC.LT.NU) GO TO 10 

0040 


NC-0 

0041 


AMP { J) — SUM/N4 

0042 


SUM-0* 0 

0043 


IF (F.GT. FM) GO TO 11 

0044 


J=J+1 

0045 

10 

F=F ♦ DF 

0046 

1 1 

DFP-N4 + DF 

0047 


SF= (DF 4- DFP)/2.0 

0048 


CALL PLOT (AHr,SF„DFP, 1 , J, 1) 

0049 

51 

CONTINUE 

0060 


RETURN 

0051 


END 
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NDATA - Number of data points in AMP array . 

Nl - Designates the first point to be plotted from the spectral 

data array by PLOT. 

N3 - Directs PLOT to skip N3 points between each plotted point 

in spectral data array* 

N4 - Number of points to be averaged when using the averaging 

feature of this routine. 

IDO - If IDO « 1 the "average 11 feature is to be used. If IDO - 3 

the "average" feature is not to be used. 

N8 “ Number of points (zeros) added by FILL. 

SUBROUTINE PRINT (NNN) 

This subroutine accepts the error cumulation array from EXT after 

a run is completed. It prints out the error table and other error data. 

Labeled Common Variables 

IE - Array containing the sum of 7 types of errors for each of 

the 300 records . 

IRCEND - The number of the last record to be processed. 

SUBROUTINE SIMP (SUM, FLL, FUL, N, NTYP , A, B, C) 

This is a Simpsons rule integration routine called by CHI, 



FUL - Upper limit. 

N - Number of points - 21. 

NTYP - Determines if routine will compute for a normal test, log- 
normal test or both. 


A - Mean value of amplitudes 

B - Log standard deviation. 


C 


Mean value of log-amplitude. 
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FORTRAN IV G LEVEL 1, HOD 4 PRINT DATE - 70122 02/95/57 


0001 


SUBROUTINE PRINT (NNN) 

0002 


COM NON/P R/I F (7 , 3Q0) 

0003 


COHMON/RAT/IRCEND 

0004 


DIMENSION IESUM (7) 

0005 


PRINT 3 

0006 

3 

FORMAT (1H1 ,50X,' DATA PROCESSING IRR EGU LAST I ES * ,/ , 1 OX , * ERROR CODES 


♦ FOLLOW* ,////, 

*10X,*NO BASE POINTS FOUND IN BASE SEARCH ********** I 1 /# 

♦ 1 0X f * NUMBER OF BASE POINTS EXCEEDS 10 ************ 2*/# 

*10X,*NUMBER OF SIGNAL POINTS EXCEEDS 10 *********** 3*/, 

♦ 1 0X , • NUMBER OF BASE POINTS INSUFFICIENT *********** 4»/, 

*10X, 'NUMBER OF SIGNAL POINTS INSUFFICIENT ********* 5*/ f 

*10X, 'NUMBER OF PASSES EXCEEDS LIMIT LI ************ B 1 /, 

♦ 1 0X NUMBER OF ERRORS IN RECORD EXCEEDS 5 ********* 7*> 


0007 


DO 555 1=1,7 

0008 

555 

IESUM {I)=0 

0009 


PRINT 1 

0010 

1 

FORMAT (10X, 'ERROR*, 11X,'1', 1 3 X , 1 2 1 , 1 3X , 1 3 ' ,13X r *4*,13X, * 5* , 1 3 X, 


**6* ,13X, *7*) 

0011 


PRINT 100 

0012 

100 

FORMAT (1 OX, * RECORD 1 ,///) 

0013 


DO 55 1= 1 , IRCEN D 

0014 


DO 55 K=l,7 

0015 

55 

IESUM (K) =IESUM (K) *IE(K,I) 

0016 


DO 2 I=1,IRCFND 

0017 


DO 12 K = 1 r 7 

0018 


IF (IE ( K , I ) ) 15, 12,15 

0019 

12 

CONTINUE 

0020 


GO TO 2 

0021 

15 

WRITE (6,5) I , (IE (K , I) r K=1,7) 

0022 

5 

FORMAT (10X,I3, 10X r I4, 1 0X, 14 , 1 0 X , 1 4 , 10X,l4, 10X,I4, 10K,I4, 10X, 14) 

0023 

2 

CONTINUE 

0024 


PRINT 20 

0025 

20 

FORMAT (10X,// r 10X,* ERROR CODE * , 1 OX , * NU MBE R OP ERRORS*) 

0026 


DO 50 K= 1 , 7 

0027 

50 

PRINT21, K, IESUM (K) 

0028 

21 

FORMAT {10X r I6, 1 4 X , 1 9) 

0029 


RETURN 

0030 


END 
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0022 

0023 

0024 

0025 

0026 

0027 

0028 

0029 

0030 

0031 

0032 

0033 
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IV G LEVEL 1, HOD 4 SIMP DATE - 70122 02/35/57 

SUBROUTINE SIMP (SUM , FLL , PUL , N, NT YP , A , B , C) 

C INTEGRAND FUNCTION REMOVE WHENCHANGING FUNCTION 

PBF (X, A, S) * (1-0/ (S*5QRT (6-28318) )) *EXP (-0.5* (( (X-A)/S) **2) ) 
FNP=N-1 

DELX= (FtIL- FLL) /F NP 

SUM=0- 0 

SUll1=0- 0 

SUM2— 0- 0 

DO 1 1=1, N 

FK=I- 1 

X = FK*DELX + FLL 

C CALL FOR INTERGRAND SUBROUTINE HERE 

C CALL FOR INTERGRAND SUBROUTINE HERE 

IF { NTYP) 101,101,110 

101 V AL=rBF (X , A , B) 

GO TO 102 

110 I F ( X) 20, 20, 100 
20 V AL = 0. 00 
GO TO 102 

100 X X= 0- 5* A LOG (X/A) 

VAL^PBF (XX ,C , 8) 

V AL -0 . 5 * V A L/X 

102 CONTINUE . 

C 

IF ( I, EQ . 1.0R.I.£Q*N) GO TO 2 
*1= MOD (1,2) 

IF ( J) 3,4,3 

3 SUfll-SUMI +VAL 
GO TO 1 

4 SUM2-SUM2+ VAL 
GO TO 1 

2 SUM=SUM ♦ V A L 
1 CONTINUE 

SDM-SUM + 2* 0*5UM 1 + 4. 0*5t)H2 

SUM=SUM*DELX/3.0 

RETURN 

END 
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SUBROUTINE FILL (I) 

This routine is called by subroutine EXT in cases where data points 
are discarded due to irregularities. FILL places zeros into the omitted 
positions. 

Argument Variables 

I - Number of data points discarded. 

Labeled Common Variables 

AMP - Array containing extracted data. 

NDATA - Number of data points in AMP. 

N8 - Number of zeros added by FILL. 


/ 
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FORTRAN 

0001 

0002 

0003 

0004 

0005 

0006 

0007 

0008 

0009 

0010 


IV G LEVEL 1 , MOO 4 FILL DATE = 70122 

SUBROUTINE FILL (T) 

COMMON/ A P/A HP (10000) f NDATA 
COMMON/NPTS8/N8 

J»U*3)/6 
N8=N8+J 
DO 1 K=1,J 
NDATA=NDATA ♦ 1 
1 AKP (NDATA) -0* 

RETURN 

END 
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ERROR CODES FOLLOW 


DATA PROCESSING I RftEGUL ART I ES 



NO BASE POINTS FOUNO IN BASE SEARCH ********** 1 

NUMBER OF BASE POINTS EXCEEDS 10 ************ 2 

NUMBER OF SIGNAL POINTS EXCEEDS 10 *********** 3 

NUMBER OF BASE POINTS INSUFFICIENT *********** A 

NUMBER OF SIGNAL POINTS INSUFFICIENT ********* 5 

NUMBER OF PASSES EXCEEDS LIMIT LI ************ 6 

NUMBER OF ERRORS IN RECORD EXCEEDS 5 ********* 7 

ERROR 1234 

RECORD 


2 0 0 0 1 

232 0 1 0 0 

240 0 0 0 l 

241 0 0 0 1 

25 2 0 0 0 I 

253 0 0 0 1 

264 0 0 0 1 

265 0 0 0 1 

276 0 0 0 1 

277 0 0 0 1 

289 0 0 0 1 

290 0 0 0 l 

300 0 0 0 1 



error COOE 
1 
2 

3 

4 

5 

6 
7 


NUMBER OF ERRORS 
0 
1 
0 

12 

0 

0 

0 


ERROR T^BLF 





TAPE NUMBER 3397 


TRACK 4 


FILE 2 


amplitude 


LOG AMPLITUDE 

COUNT 


CUMULATIVE PROBABILITY 

0 • 235000E 

03 


-0.347886E 00 

0* tOOOOOE 

01 

0.120077E-03 

0.245000E 

03 


-0.327050E 00 

0-200000E 

01 

0* 360230 E-03 

0 • 255000E 

03 


-0.307047E 00 

0*2000006 

01 

0.600384E-03 

0*265 000 E 

03 


-0.287Q14E 00 

0*0 


0.600384E-Q3 

0 *2750006 

03 


-0.269293E 00 

0*0 


0*6003846—03 

0 * 285000 E 

03 


-0*2514346 00 

0 • 200000E 

01 

0* 840538E— 03 

0 * 295000 E 

03 


-0.234191E 00 

0*1000006 

01 

0*96061 5E—03 

0* 3G5000E 

03 


— 0 *2175238 00 

0* 300000E 

Oi 

0* 132G85E— “02 

0*31 5000E 

03 


-0.20L393E 00 

0.600000E 

01 

0* 204I31E— 02 

0*325000E 

03 


— 0 ♦ 1 8 5766E 00 

0*8000006 

01 

0*3001926-02 

0 • 335000E 

03 


-0.L70614E 00 

0.4000006 

01 

0.348223E— 02 

0*3450006 

03 


— 0 .1559076 00 

0.170000E 

02 

0.552353E-02 

0.3550G0E 

03 


-0.14L620E 00 

0 * 190000E 

02 

0* 780499E— 02 

0* 365000E 

03 


-0.127730E 00 

0*2100006 

02 

0* 103266E— 01 

0*3750006 

03 


-0.114216E 00 

0*390000 6 

02 

0*1500966-01 

0* 385000 E 

03 


-0.10I057E 00 

0.660000E 

02 

0.229347E-01 

0 * 395000 E 

03 


— 0*33235 9E -01 

0*7200006 

02 

0*3158026—01 

0 • 405000 E 

03 


-0. 7573536-01 

0 * 127000E 

03 

0*4683006—01 

0.4150G0E 

03 


-0*6353966-01 

0 .1780006 

03 

0* 682036E— Oi 
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APPENDIX B 

This program generates a set of N random numbers having a log- 
normal distribution and a pre-selected mean and standard deviation. The 
program is in the form of a FORTRAN IV subroutine. 

Theory : By definition a log-normal random deviate is one whose logarithms 

are normal random deviates. Thus if (X^) is a set of log-normal random 
numbers then there must exist a set of normal random numbers (y^) related 
to the X^ by 

y^ = In X^ B1 

Equation Bl may be generalized by the addition of appropriate scaling 
factors; i.e., we may let 


y. = a In X, + b B2 

Now by choosing the mean and variance of the (y^) and the values of the 
scale factors a and b it is possible to generate a set of (X^) having any 
desired mean and variance from a set of normal deviates (y ) . Solving 
B2 for we have 

y - b 

X. - exp (— ) B3 

1 a 

Since we wish to specify only two parameters, viz,, the mean and standard 

deviation of the (X^) it seems reasonable to assume that we will need 

only two parameters in equation B3. We therefore let a = 1 and take the 

mean of the (y.) to be zero. B3 then becomes 
x 

X ± = exp (-b) exp (y ) 


B4 



101 


taking the average of both sides of equation B4 we have 

X * exp (-b) exp(y^) B5 

and also taking the second moment of (X^) about zero 

X 2 = exp (-b) exp (2y^) B6 

the averages of the exponential functions in equation B5 and B6 can 
be evaluated easily 


exp (ny i ) 


(2irt 2 ) 


- 1/2 


exp(ny) • exp(y /2a)dy 


B7 


-x 


Combining equation B5, B6 and B7 we obtain expressions which may be 
solved for the scale factor b and the required standard deviation of 
the (y i ) 


o 2 = In (p/X 2 ) 


B8 


and 


exp (-b) 


- Xexp(-o /2) 


B9 


where y is the second moment of the (X^) about zero. 

Program: The log-normal generator makes use of the normal random number 

generator included in the IBM Scientific Subroutine Package for the 360 
computer. This routine (GAUSS) generates normal random numbers with any 
required mean and standard deviation. Coding for the program is shown 
in the accompanying listing. The argument list is as follows: 

AVE - The required mean. 


VAR - 


The required standard deviation. 
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Y - A vector of log-normal random numbers returned by the 

subroutine. Y is dimensioned by the calling program. 

N - The number of random numbers to be generated. 

IX - A "seed" required by GAUSS. IX must be a 5 digit odd integer. 

Statements 3 to 6 compute the required standard deviation for the 
Gaus sian-random numbers and the proper scaling factor. Statements 7 to 

9 call GAUSS compute a log-normal random number from equation B5. 

Fortran List for Log-Normal Generator 

1 SUBROUTINE LOGN (AVE , VAR, Y, N, IX) 

2 DIMENSION Y(l) 

3 VAR = VAR* *2 + AVE **2 

4 SIG = ALOG ( VAR/ AVE**2 ) 

5 Z BAR = EXP (SIG/2. 0) 

6 SIG = SQRT(SIG) 

7 DO 1 I = 1, N 

8 CALL GAUSS (IX, SIG, 0.0, X) 

9 1 V(I) = (AVE/Z BAR) *EXP (X) 

10 RETURN 


11 


END 
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Appendix 4 


REQUIREMENTS STUDY FOR THE AVLOC EXPERIMENT 

This Report describes the various analyses which have been made 
in connection with the MSEC High Altitude Aircraft Test. These analyses 
include an investigation of the effects of the aircraft navigational 
errors on experimental accuracy, a survey of the engineering measurements 
to be conducted on this experiment and a survey of the meterological 
support requirements for the project. In addition a brief description of 
the High Altitude Aircraft Test program is given as background for the 
reader who may not be familiar with this project. Volume Two contains 
a complete report of both the experimental and theoretical investigations. 
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DESCRIPTION OF THE MSFC HIGH ALTITUDE AIRCRAFT 
TEST FOR VISIBLE LASER COMMUNICATIONS 


Introduction 

The analysis reported in the following sections have been performed 
in support of the High Altitude Aircraft Test of Visible Laser Communi- 
cations being conducted by MSFC personnel during 1971-1972. In order 
to orient the reader as to the relation of this work to the overall 
MSFC Aircraft Test program a brief description of that experiment 
will be given. A more comprehensive discussion of the experiment can 
be found in the High Altitude Aircraft Test project plan [2]. 

Program Objectives 

The High Altitude Aircraft Test for Visible Laser Optical Communi- 
cations (hereafter referred to as HAAT) is intended to perform three 
functions, viz., to collect scientific data on the propagation of 
visible wavelength radiation through the atmosphere, to provide engi- 
neering data needed for the evaluation of the techniques of optical 
communications and for the design of future systems, and to demonstrate 
the feasibility of optical techniques for communications between air- 
craft and ground or satellite and ground. These techniques include 
the ability of the satellite (or aircraft) to acquire the ground 
station and vice versa, the ability to accurately track the moving 
satellite and the ability to transmit high data rates with low error 
and high reliability. 
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Experimental System 

The experimental system consists of an airborne optical communi- 
cations terminal carried aboard a modified RB-57 aircraft and a ground 
station located on Redstone Arsenal, Alabama. Both the airborne and 
ground terminals are equipped with transmitting lasers so that two way 
communications may be established and both uplink and downlink 
propagation may be studied. In addition to communication equipment 
each terminal contains an optical tracking system which controls its 
pointing during the experiment. Thus each end of the system coop- 
eratively tracks upon the other end throughout the course of the 
experiment. The ground station is equipped with a separate laser 
radar (designated the Ground Based Acquisition Aid or GBAA) for 
initial acquisition of the aircraft and for providing exact range and 
zenith angle information throughout the experiment. Both the airborne 
and ground terminals are fully instrumented to record all pertinent 
data relating both to system operation and to atmospheric effects on 
optical propagation. 

o 

The downlink channel will consist of a 6328 A He-Ne laser beam 
pulse code modulated at 30 M. bits /sec. Downlink transmission will 
consist of either a real time television picture, a pseudo random 
code word for bit error rate measurements, or telemetry of data 
being collected aboard the aircraft. The uplink channel will consist 

o 

of an amplitude modulated 4880 A Argon laser beam which will be used 
for uplink scintillation measurement and to transmit commands to the 
airborne package. 



The ground .terminal receiving and transmitting optics consists of 
a 2 4- inch Cassegrainian telescope with the experimental package located 
at its Coude focus. The incoming 6328 A beam from the airborne He-Ne 
laser passes through a wavelength selective beam splitter to separate 
it from the 4880 & uplink beam. It is then directed by means of a 
second beam splitter to two detectors. The first detector is a quadrant 
photomultiplier which provides fine pointing signals. These signals, 
along with the output of the angle encoders on the telescope mount, are 
fed into a SCC 4700 computer which performs the necessary coordinates 
transformations and generates the signals to drive the polar-equitorial 
telescope mount. In addition, both the angle encoder and the quadrant 
photomultiplier outputs are recorded to give angle of arrival fluctu- 
ation data. 

The second part of the incoming beam is directed to two photo- 
multiplier tubes. The output of one of these detectors is recorded to 
provide the wide bandwidth downlink communication, channel output, 
scintillation data, and bit error rate data. 

The uplink beam is generated by an Argon laser which passes 
through a modulator and a beam steerer and out through the telescope. 

The beam steerer is driven by the computer so as to direct the uplink 
beam directly back along the incoming beam. The uplink beam is 
modulated with a 10.7 MHz sine wave and also with audio tones used to 
control the operation of the aircraft terminal. 

In addition to the main tracking and communications system, the 
ground terminal contains a laser radar for acquisition and ranging 
(the GBAA) . The acquisition radar is mounted on the main telescope 
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tube and bore-sighted with it. It consists of pulsed Argon laser, 
beam steering optics and receiving electronics. During acquisition 
phase, the beam is digitally scanned in a raster pattern, until a 
return from comer reflectors mounted on the aircraft is detected. At 
this time the GBAA enters a limited scan mode and the tracking and 
pointing function is passed to the detector in the primary receiver 
system. The GBAA continues to track the aircraft, however, and to 
provide range information. In the event that the primary tracker should 
lose the target, the acquisition radar will automatically re-enter the 
acquisition mode and reacquire the target. During this time the tele- 
scope drive will enter a coast mode so as to continue to point at the 
assumed position of the aircraft. 

The airborne optical system will be mounted in a fixed position 
aboard the RE-57 aircraft and will view the ground in a steerable 
mirror which will be servocontrolled to point the outgoing beam in the 
proper direction. During acquisition the steerable mirror is pointed 
in the general direction of the ground station manually. When the 
ground based acquisition radar illuminates the aircraft, an acquisition 
sensor detects the upcoming beam and causes the system to enter a 
track mode. The incoming beam passes through a dichroic beam splitter 
which isolates it from the outgoing beam to an image dissector which 
provides tracking information and a photomultiplier which detects up- 
coming commands and measures scintillation. The transmitter section 
of the airborne terminal consists of a He-Ne laser which is super- 
imposed on the incoming beam by means of the dichroic beam splitter. 

The outgoing beam is pulse code modulated at 30 megabits /sec with 
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either a pseudo random word for bit error rate measurements, with a 
video signal generated by a television camera aimed at the ground, or 
with telemetry. 

Experiment Plan 

The experiment plan calls for four missions at various times during 
the year. Each mission will consist of four fly overs of about three 
hours duration. 

The aircraft will approach the Madkin Mountain ground station on 
a straight line path at an altitude of about 50,000 feet. Once 
acquired by the GBAA radar the aircraft will assume a circular flight 
path centered on the ground terminal thus maintaining a constant range 
from the ground station and constant zenith angle. Initially an 
altitude of 70,000 feet and a path diameter of 10 miles is planned. 

This corresponds to a start range of 74,800 feet and a zenith angle 
of 20.6 degrees, the minimum zenith angle obtainable with a circular 
path. By varying the altitude of the aircraft and the diameter of 
it f s flight path the range and zenith angle can be varied at will. 

The scientific measurements to be made during each fly-over are 
outlined in Table I. Detailed discussions of these experiments can 
be found in the referenced measurement program document [2]. 



Table I . Measurements Outline 


Quantity Measured 

1, Scintillation 
(downlink) 


2, Scintillation 
(uplink) 


3. Angle of Arrival 
Fluctuations 
(downlink) 


4. Angle of Arrival 
Fluctuations 
(uplink) 


5. Bit Error Rate 


Parameters 


Analysis to Yield 


Receiving Aperture 
Range 

Zenith Angle 


Log amplitude variance 
Probability density function 
Verify theoretical predic- 
tions concerning zenith 
angle , and range dependence . 
Aperture averaging effects: 

a. Reduction of signal 
variance 

b. Change of probability 
density function. 

Difference in up and down 
link. 


Receiving Aperture 
Range 

Transmitter Aperture 


Log amplitude variance 
Probability density function 
Verify theoretical predic- 
tion concerning zenith angle 
and range dependence. 

Verify no uplink aperture 
averaging 

Effect of transmitting 
aperture size. 

Differences in uplink and 

H nun 1 *f nV . 


Range 

Zenith Angle 
Receiving Aperture 


Range 

Zenith Angle 


Range 

Zenith Angle 
Beam Divergence 
Transmitter Power 


Variance 

Probability density function 
Aperture averaging 
Dependence on range and 
zenith angle. 

Differences between uplink 
and downlink. 


Variance 

Probability density function 
Dependence on range and 
zenith angle. 

Differences between uplink 
and downlink. 


Verify theoretical predic- 
tion of BER dependence on 
system noise and irradiance 
fluctuation variance. 
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Table I. Measurements Outline (Cont’d) 


Quantity Measured 

Parameters 

Analysis to Yield 

6. Atmospheric 
Transmittance 

None 

Atmospheric transmittance 
at optical wavelengths. 

7 . Engineering 
Measurements 


Determine ability to 
acquire and track. 
Evaluate system 
performance. 



ANALYSIS OF THE ERRORS IN THE MSFC OPTICAL COMMUNICATIONS 
EXPERIMENT DUE TO FLIGHT PATH INACCURACIES 


Introduction 

Marshall Space Flight Center's High Altitude Aircraft Test of 
Visible Laser Communications (HAAT) is expected to provide scientific 
data concerning atmospheric turbulence and its effect on the propagation 
of laser beams over near vertical paths. If these data are to be 
meaningful, it is necessary that certain parameters such as aircraft 
range and zenith angle be held constant. Clearly no aircraft can fly 
a precise path, so some variation in these parameters must be accepted. 
It is the purpose of this report to examine the errors introduced into 
the experimental data as the result of deviation of the aircraft from 
its prescribed flight path and, based on this analysis, to suggest the 
maximum deviations which can be tolerated. 

General Considerations 

A description of the HAAT experiment is found in the HAAT project 
plan [2], During the test the aircraft will fly a circular pattern 
over the ground station, ideally maintaining a constant range and zenith 
angle. The ground station will track the aircraft by means of a ground 
based acquisition aid (GBAA) also used for initial contact. The GBAA 
is an optical laser radar with ranging facilities and will provide 
continuous range and angle information. From the GBAA data the air- 
craft’s position relative to the ground station will be accurately 
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known, a posteriori. It is therefore not a requirement that the air- 
craft fly precisely along a predetermined path but only that it remain 
upon a circle reasonably close to the desired pattern. It is the 
variations of the aircraft path from this somewhat arbitrary circle 
which are important. The GBAA data will also allow the elimination 
of any data segment for which the variation in the range or zenith 
angle are excessive. 

If the aircraft were to fly a perfectly level, circular path then 
its range (r) and zenith angle (0) would be constant. We anticipate 
that the aircraft 1 s actual path will be some smooth curve which approxi- 
mates the prescribed circular pattern. We will therefore assume that 
0 and r will be slowly varying functions of time. 

The quantities to be measured (such as scintillation) are statisti- 
cal quantities whose mean and variance depend upon r and 0 as well as 
upon the statistics of the atmospheric turbulence. Since the atmosphere 
itself is a non-statlonary system, the time variation of r and 9 will 
simply serve to increase the variation in the statistical parameter of 
the measured quantity. It therefore seems reasonable to require than 
the variation in the statistical estimates of these parameters due to 
aircraft motion be small compared to the variations introduced by such 
unavoidable effects as atmospheric non-st ationarity. 

Table I in Appendix B of the HAAT program plan [3], which has been 
reproduced as Table I of this report, describes seven measurements 
which will be made. Eliminating from consideration the Engineering 
measurements and neglecting any difference in uplink and downlink 
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measurements, these may be reduced to four basic types of measurement 
viz. 

1. Scintillation 

2. Angle of Arrival Fluctuations 

3. Bit Error Rate Measurements 

4. Transmittance. 

Each of these measurements will be considered. 


Scintillation 

The variance of the log-amplitude c ^ for a spherical wave is given 
by [4] 


a 


% 



r 2 (\,\ r ( 2 ~ 5 ) S 1 
C N (h) 1 


5/6 


ds 


(i) 


where C 2 (h) is the index of refraction structure constant at an 
N 

altitude h, z is the slant range and ds is an element of length along 

the propagation path* The actual situation which we wish to consider 

is that of a collimated, truncated gaussian beam. Since, however, we 

are only concerned with the order of magnitude of the variation of 

2 

on range and zenith angle, and since the form of the function C N (h) is 

not well known, we will assume that the relation of equation (1) 

describes the gaussian beam accurately enough for our purposes. 

2 

We will assume a form for C N (h) viz. 


2 2 -1/3 ' h/h 0 

V h) = V h e 


( 2 ) 


where h Q is a scale height usually taken to be 3200 m. Recent work [5] 

indicates that this value may be too large, however since for larger 

h , a 0 is more sensitive to changes in altitude and zenith angle we 
u * 

will assume this value as a worst case. 
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For a zenith angle 0, h is related to 3 by 


h = s cos 9 


( 3 ) 


then 


2 _ 2 -1/3 . 

a JL = C 0 COS 6 


2 

l 


{ s cos 6 \ 

S ~ 1/3 e h ° [ ] 5/6 ds (4) 


Now we let u = s/z, then 


a 


2 

S, 



U 1/2 (1-U) 5/6 


-u z cos 


0/h 


du 


(5) 


But z cos 0 is just the aircraft altitude h , therefore 
^ m 



cos 


- 11/6 3/2 


m 


I 


„ 1 / 2 a - U ) 5/6 



du 


( 6 ) 


Now consider an aircraft flying in a circular path of radius R at an 

altitude h . As can be seen in Figure 1, if the radius of the circle 
m 

changes by an amount AR then h^ remains constant and 0 changes by an 
amount A0 given by 


A0 


AX _ AR cos 9 
s s 


h AR 
m 


( 7 ) 


then the fractional change in the log amplitude variance is 
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Ac 


2 

J. 



1 


2 


G 



0 


_ r+ ii c : 

" l+ 6 C 0 


-17/6 

cos 0 


sin 


G h 


3/2 


m 


1 <H ) ] 

m 


A0 


/V cos e I(h n ), <»> 

Where I(h ) denotes the definite integral in Equation 6* Equation 8 
m 

reduces to 

2 

E = m ■— tan 6AG (9) 

AR 2 6 

°£ 


or 


11 RAR 

■ " v 


( 10 ) 


Equation 10 allows us to estimate the variation in for a 

change in horizontal position of the aircraft. A similar equation may 
be derived for a change in altitude AH. As can be seen from Figure 2, 
a change in aircraft altitude changes both h m and 6. The change in 6 
is given by 


then 


A0 = 


AX 


A H sin 6 
m 


R AH 


m 


(ID 


Ao, 


Bo, 


30 


A0 + 


dh m 
m 


( 12 ) 


i 
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Ac, 


■4- tan e A0 + 4 h 1 Ah + i" 1 ^ ) 

6 2mm m 


91 ( h J 

nt 

3h 

m 


Ah 


m 


(13) 


Aa n 


H 


11 R_ 

6 s 2 


2 AH 


m 


H 


m 


* I 


AH 


m 


H 


m 


1_ 

I 


m 


31 (h ) 
m m 

9h 

m 


AH 


m 


(14) 


The last term in Equation 14 has been evaluated numerically on the 
IBM 360/50 computer using a double precision Simpson's rule integration. 
The values of this term are 

- 2.45% per 1000 ft. at 42,000 feet 

- 2.15% per 1000 ft. at 52,500 feet 

- 1.90% per 1000 ft. at 63,000 feet. 

2 

Equation 10 and 14 specify the approximate errors in a due to 
variation in AR and AH, These have been evaluated for the altitudes 
indicated above and for horizontal distances of 5 and 10 miles. The 
results are shown in Table II. The values for AE^ are seen to be 
quite small and may have either sign. This behavior may be understood 
by considering that for a positive AH, corresponding to an upward motion 
of the aircraft, the log amplitude variance is increased due to an 
increase in the propagation path length. 


Table II 


Altitude 

Range 

Slant Range 
(1000 's ft.) 

(1000’ s 

ft) (Miles) 

42.0 

5 

49.2 

52.5 

5 

58.8 

63.0 

5 

68.3 

42.0 

10 

67.2 

52.5 

10 

74.5 

63.0 

10 

82.1 


a e r 

a e h 

| AE j + 1 AH | 

%/1000 ft. 

2/1000 ft. 

%/1000 ft 

2.0% 

-.11% 

2.1% 

1.4% 

+ .005% 

1.4% 

1.04% 

+ .05% 

1.1% 

2.15% 

.85% 

3.0% 

1.75% 

1.11% 

2.9% 

1.44% 

1.17% 

2.6% 
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2 

This effect is relatively slight; however, since C N at the alti- 
tudes we are considering is very small. Furthermore, it is compensated 
for by the decrease in zenith angle which tends to swing the propa- 
gation path up out of the more turbulent lower atmosphere. The variation 
in the horizontal distance has a larger effect since a positive AR (the 
aircraft moves outward from the ground station) both lengthens the 
propagation path and lowers it into a more turbulent region. 

Angle of Arrival Fluctuations 

The variance of the apparent angle of arrival can be shown to 
depend on the five-thirds power of the correlation length ^[6.7], i.e., 

<> 2 > cc l/r o 5/3 (15) 

For a spherical wave propagating from the ground to an aircraft at a 
slant range z. r^ may be expressed as 



or 


■)Z cos 0/h 


l/r 5/3 * 
o 


,, v -1/3 -u , U h o 

(h o u) e ( Z COS 0 


5/3 


■) 


h du 

<— rr) a?) 

cos □ 


where we have let u — s cos 0/h^« We note that z cos 6 is the aircraft 

altitude h , therefore 
m 


<«*> 


oc 


7/3 


h 5/3 cos 0 
m 


h /h 


m 


-u 4/3 . 
e u du 


(18) 
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The integral in Equation 18 is the incomplete garoma function r(h/h Q ,7/3) 

which, at the altitudes that we are considering, is a slowly varying 

function of h . We will therefore take it to be a constant and write 
m 

/Ae 2 ^ of h _5/3 cos -1 9 (19) 

For changes in horizontal range 


E,„ - A ^ - Tan 6 46 (20) 

“ <46 2 > 


or 



( 21 ) 


By comparing this expression with Equation 10, we see that the 
angle of arrival measurement errors are smaller than the scintillation 
measurement errors by a factor of 6/11 for changes in R. Hence the 
requirement on the accuracy of R is not determined by the angle of 
arrival measurements . 

For changes in the aircraft altitude we may replace cos 0 by 

h /s hence 
m 

< 46 2 > cc h - ?/3 . - h - 7/S t/r 2 + h 2 1 (22) 

' ^ m m y m 

then 


E AU = (~ i H 1 + H s" 2 )AH 
AH 3 m m m 


(23) 


For a five-mile radius E„ varies from 3.8% at 42,000 feet to 

n 


2.5% at 63,000 feet. 



18 


For downlink angle of arrival , the integral in Equation 16 is 
replaced by 



C>) <«, 


5/3 


ds 


No calculation for downlink angle of arrival fluctuation has been 


made. 


Bit Error Rate Measurements 

The probability of error for a binary coded channel in the presence 

of log normal scintillation has been given by Fried [8] as 

, , - - 1/2 
P E = [2tt <8B JCj +p Xj)] 

x [exp {- -j (48^0^ Xj, +S^ X T > } ] (24) 

2 

where B is the signal to noise ratio in the receiver and a is the 
log amplitude variance of the scintillation. The parameter X^, is 
related to o^ and B by 

a* = -(in X t )/(4 + 8B 2 Xj) (25) 

Equation 24 was derived under the assumption that the receiver noise 
is white gaussian noise. This assumption may not be justified for an 
optical receiver where the principle noise source is the Shott noise 
in the photodetector. A more precise treatment of detection in an 
optical communications system has been given by Hoversten [9]. For 
the purpose of estimating the dependence of on the small changes in 
range the gaussian noise model should be adequate, however, so that 
Equation 24 may be used. 
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Inspection of Equation 24 reveals that the probability of error, 

P , depends only on the signal to noise ratio 6 and the log amplitude 
E 

n 

variance of the scintillation . Provided that we may assume that the 
receiver noise remains constant then 6 is proportional to the received 
power which in turn varies as the square of the slant range. Hence, 

M = ? — ( 26 > 

6 S 

For a 1000 foot change in slant range, a 5 mile radius path and altitudes 
of 42,000 and 63,000 feet this error amounts to 4.76% and 2.44%, 
respectively . 

The changes in a. have been discussed in a preceeding section. 

Xj 

It therefore remains only to demonstrate how these changes will effect 
P E* 

Equation 24 may be simplified by noting that for a small change in 
6 or o 2 the first factor will be slowly varying compared to the 
exponential. We may therefore treat it as constant. Eliminating 
we ob tain 

, 46 4 InL ■ L 2 

P =Kexp[-i( -+B V ] (27) 

E z 4 + 86 Xj, 

or 

Pg = K exp [-fCXj.g)] (28) 

then 

_ i_ 9 p e _ a f 

P E 3 *T * ' 3 *T 


( 29 ) 
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and 



9f 

9e 


(30) 


From Equation 25 it is seen that may take on values from a 
very small positive number (deep scintillation) to +1 (no scintillation) 
while $ can take on any positive value. For the case of fairly large 
signal to noise ratio and small scintillation f(B,X r[ ,) will reduce to 

f (B ,X T ) = j B 2 (Xt + J in Xy) (31) 


then 


Then 


3 f _ 1 /i , 1 \ a/ ,5. 

9 Xj " 2 6 (1 + 2X t ) * / 4 6 


|| - (x T+ |£nx r )- 



9 Xj. 



!is 
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BX^B 



if 5^ - B V e 


(32) 


(33) 


(34) . 


(35) 


(36) 


Now for the assumed conditions, i.e., B large and approximately 
equal to unity the first term is small so that 



sx t ab 


(37) 
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From this analysis it is clear that even these with the simplifying 
assumptions the bit error rate may be a strong function of the signal 
to noise ratio, particularly for the case of large It may therefore 
be necessary to limit the experimental conditions to those producing 
large bit error rates. 

A more exact analysis could be carried out by avoiding the 
simplifying assumptions which we have made and differentiating 
equation directly * However, due to the uncertainty as to the accuracy 
of the gaussian noise model it was not felt that such an analysis 
would be worthwhile. 


Atmospheric Attenuation 

The atmospheric transmission is given by [10] 

Tr(q S) 2 I 

t a P m 


(38) 


where S is the slant range to the aircraft; a the laser beam divergence 
angle, the mean transmitted power and I the mean received power. 

In the preceeding section it was shown that for an altitude of 63,000 
feet, I would change by about 2-44% per 1000 feet change in slant range 
For attenuation measurements this variation is averaged so that small 
changes in the slant range have little effect on the observed value 
of t provided that the average value of the slant range is accurately 

A 

known. This information should be available from the GBAA. We there- 
fore conclude that for changes in aircraft position of the order of a 
thousand feet will introduce errors of only a fraction of a percent 
into the measured value of x , . This is far less than the errors which 
may be expected in the measurement of I and P^ due to inaccuracies 
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in the calibration of the detectors. 

The atmospheric extinction constant K is related to by 

K » - nr 2.n T (39) 

M A 


Where M is the air mass. M is proportional to 



(40) 


Where S is the slant range and h is height above the ground. For the 
altitudes we are considering this integral is constant to within about 
1 / 2 %. 


Conclusions 

From the preceeding analysis , it is clear that the quantities to 
be measured on the HAAT experiment are relatively insensitive to small 
changes in either the altitude of the aircraft or its horizontal dis- 
tance from the ground station. Of the two the latter is the more 
critical. In general, changes in aircraft position of 1000 feet will 
produce variations of only 1 or 2% in the measured statistical param- 
eters. A typical data segment will probably be 1 or 2 minutes long. 
Experimental results indicate that the measured values of log ampli- 
tude variance taken a few minutes apart may differ by a considerable 
amount. It is therefore not unreasonable to expect the log amplitude 
variance to change by a few percent during the course of a measurement. 
If this is the case, then an additional variation of a few percent 
due to flight path inaccuracies will not be objectionable. 
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A possible exception is the bit error rate measurements. At 
very small bit error rates, the error rate becomes a very sensitive 
function of signal to noise ratio and therefore a sensitive function 
of received power. This, combined with the fact that a long time 
is required to measure a very small bit error rate, may make it 
difficult to compare the measured BER with theory. 

Recommendations 

It is recommended that the initial operational plan specify 
navigation tolerances during data acquisition of +500 feet in range 
from the ground station and +500 feet in altitude. These tolerances 
are to be held for one complete orbit around the ground station. 

If it appears likely, either from discussion with Air Force 
personnel or from experience on early flights in the program, that 
these tolerances Impose excessive restraints upon the operation of 
the aircraft, then they may be relaxed by as much as a factor of 2 or 
3 without jeopardizing the validity of the experiment. 



METEROLOGICAL DATA REQUIREMENTS 


Introduction 

The High Altitude Aircraft Test for Visible Laser Optical Communi- 
cations requires that meterological data be collected prior to and 
during each test flight so that the optical measurements may be corre- 
lated with such gross atmospheric parameters as temperature, barometric 
pressure, relative humidity, and wind profile. In addition meteorological 
data will be useful in revealing unusual atmospheric conditions, such 
as a temperature inversion, which might lead to anomalous experimental 
results. A second, and equally important, requirement for meteoro 
logical data is to assist is scheduling the test at the most advantageous 
times from the point of view of favorable weather conditions. Also the 
test conductor must be provided with current and accurate forecast so 
that in the event of deteriorating weather conditions he will have 
the best possible information on which to base a decision to cancel or 
postpone the test flight. 

Several sources of meteorological data are available. These include 
instrumentation which will be located at the Madkin Mountain ground 
terminal, the MSFC Atmospheric research station facilities located 
approximately three miles southwest of the ground terminal, the FAA 
Aircraft Advisory westher facilities located at the Huntsville - Decatur 
Airport five miles to the west, and several other weather stations on 
Redstone Arsenal operated for other purposes. In addition use will 
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be made of the normal NOAA wide area weather information services. 

The MSFC Atmospheric Research Station will be the principal source 
of meteorological measurements and will assume primary responsibility 
for the collection and evaluation of weather information from all 
sources and for providing weather forecasts to the test conductor. 

We have, at the direction of the contracting officers repre- 
sentative (COR) , undertaken to review the meteorological support 
requirements for the High Altitude Aircraft Program, A preliminary 
statement of the meteorological data to be collected in connection 
with these test is given in the HAAT program plan [11]. We have care- 
fully reviewed the meteorological measurements specified in that docu- 
ment to determine their value in analyzing the optical propagation 
data which will be obtained and to insure that no additional metero- 
logical data will be required. Several conferences have been held 
between the project director and MSFC Atmospheric Research Station 
personnel. On the basis of these conferences and our review, specific 
recommendations for meteorological data requirements have been formulated. 
We have also considered the question of what, if any, specific criteria 
can be established for the minimum weather conditions under which a 
flight test will be conducted. 

Application to Measurements Program 

Meteorological data will be used in the analysis of both the opti- 
cal propagation experiments and the engineering experiments performed 
during the high sltitude aircraft test. Areas in which meteorological 
data will be employed include the following. 
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a) General Weather Conditions » For the optical propagation data 
to be meaningful it is necessary that the general weather conditions 
under which it was collected be carefully and quantitatively documented. 
The same is true, perhaps even to a greater degree, for the analysis 

of the communication and tracking systems performance. For this 
reason careful meteorological observations must be made throughout 
each flight. Special importance must be attached to the identification 
of any unusual conditions, such as temperature inversions or high, thin 
cloud cover, which might not be readily apparent to a casual observer. 

b) Wind Profile . The frequency spectrum of the optical scintil- 
lation depends upon the component of wind velocity transverse to the 
optical path. Wind profile measurements will therefore be directly 
applicable to the analysis of the scintillation data. Fried [12] has 
developed expressions for the scintillation spectrum as a function of 
wind speed assuming that the wind speed is constant along the optical 
path. To our knowledge no analysis has been performed for the case of 
wind speed which varies along the optical path. Therefore additional 
theoretical work may be required in this area in order to properly 
interpret the results of the aircraft test. 

c) Estimation of Structure Constant Profile . In order to ade- 
quately compare the experimental results to theory it would be desirable 

2 

that the atmospheric structure constant, be measured as a function 

2 

of altitude during each test. Methods of determining from the 
temperature structure constant have been proposed by Ochs [13,14]. It 
does not appear, however, that the instrumentation to make these 
measurements capable of being flown on a weather balloon will be 
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available in time for use on this project. It will therefore be 

2 

necessary to assume a standard profile for and to attempt to infer 

from conventional meteorological data the validity of this assumption. 

Information as such factors as temperature lapse rates will be useful 

2 

in detecting conditions under which the profile may deviate from 
normal. 

Test Scheduling 

One of the most important applications of meteorological infor- 
mation will be it's use in scheduling the test flights at times when 
the weather conditions are likely to be the most favorable and in 
deciding when conditions warrent the postponement or cancellation of 
a test. The problem of scheduling will consist of three phases, viz. 
(1) Advanced scheduling of flight periods, (2) selection of specific 
flight times and (3) the decision to delay or abort a test in the 
event of deteriorating weather conditions. 

We recommend that the following procedure and criteria be adopted 
for scheduling of flight test from a weather standpoint and for 
cancelling of a mission once it has begun. 

A. Advanced Scheduling 

The test will be carried out in sequencies of four flights. The 
operational periods for each sequence will be determined as far in 
advance as possible. In order to schedule the operational periods 
at times when most favorable weather conditions are likely to prevail 
it will be necessary to compile statistical data concerning the 
typical weather patterns at various times of the year in the Huntsville 
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area. Thi 9 data will include such factors as the average number of 
clear days per month, the typical monthly rainfall, etc. Based on 
statistics available for the past several years preferred operational 
periods will be established. 

B. Selection of Specific Flight Times 

One week before the beginning of each operational period tenta- 
tive flight schedules will be established- The specific dates and 
times chosen will be based primarily upon a prediction of a low percent 
cloud cover and the absence of fog* haze or other conditions which 
might adversely effect the operation of the experimental system* The 
schedules must be considered tentative at this point since the 
reliability of one-week forecast is not very great. 

Flight schedules will be reviewed at 5 days, 3 days, 48 hours and 
24 hours before the beginning of each test. Rescheduling will be 
possible at each of these points if revised forecasts indicate an 
increased probability of fog, haze or excessive cloud cover. 

C. Cancellation of Test 

For the 24 hour period prior to the flight revised forecasts 
will be made every 6 hours and then hourly for the 6 hour period 
immediately preceeding the flight. 

D. Criteria for Postponement of Test 

It is recommended that the following criteria be used in scheduling 

and delaying a test flight - 

1. Predicted cloud cover greater than 10 percent. 

2. Any prediction of fog, haze, or cirrus clouds which 
might reduce atmospheric transmission below 70%. 
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Data Requirements 

On the basis of our survey of the meteorological data required to 
support the High Altitude Aircraft Test and as the result of several 
conferences held with the MSFC meteorologist we recommend that the 
following meteorological measurements be made. 

A. Madkin Mountain Ground Station 

A meteorological shelter will be installed near the Ground Terminal 
on Madkin Mountain. Instrumentation will consist of a hydrothermograph 
to provide a continuous recording of relative humidity and temperature 
and a recording anamometer to measure wind speed and direction. Measure- 
ments at this site are required from 24 hours prior to the beginning of 
each test flight until 24 hours after the conclusion of the test. 

B. At MSFC Atmospheric Station 

1. Hourly observations from 6 hours prior to the beginning of each 
flight test until 1 hour prior to the beginning of the test, and then 
every 1/2 hour until the conclusion of the test. Observations will 

consist of: 

a) Temperature 

b) Barometric pressure 

c) Relative Humidity 

d) Wind Speed and Direction 

e) Clouds, percent cover, type and approximate altitude 

f) Surface visibility 

2. Radiosonde Releases to provide temperature, relative humidity, 
wind speed and wind direction as a function of altitude from the surface 
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to 110,000 feet at 1000 foot intervals. Relative humidity measurement 
will not be made above 40,000 feet due to its low value above this 
level. It is expected that all of the sondes may not reach the speci- 
fied altitude. This is not critical, however, since data above the 
operating altitude of the test aircraft is of only minor interest. 

Two radiosonde releases are required for each flight test. The 
first will be made enough in advance of the arrival of the aircraft 
over station so that a complete profile will be available at the 
beginning of the test. The second release will be made near the 
mid-point of the test. In the event that the aircraft remains on 
station for the maximum time (4 hours) a third release near the end 
of the test would be very desirable. 

C. Wide Area Weather Information 

To aid in evaluating the local weather data and in predicting 
conditions at the time of each test wide area weather information is 
required. Meteorological FAX charts will therefore be supplied for 
surface, 500 mb., 700 mb. and 850 mb. every six hours of 48 hours 
preceeding and 24 hours following each flight. 



ENGINEERING MEASUREMENTS PROGRAM 


Introduction 

The High Altitude Aircraft Test program plan contains, as Appendix B, 
a detailed measurements program for the scientific experiments which 
will be conducted. Engineering measurements are briefly outlined in 
Section III.C.5 (page 22) of the subject document. It was desired by 
the cognizant MSEC personnel that the engineering aspects of the HAAT 
program be documented in detail comparable to the scientific measurement 
program. We have, therefore, prepared at the request of the contracting 
officers representative this engineering measurements plan which describes 
the engineering and system analysis aspects of the Aircraft Test program. 

Engineering Measurements 

The measurements described in Sections III.C.l through III.C.4 of 
the Measurements Program for the High Altitude Aircraft Test of Visible 
Laser Optical Communications are intended to yield fundamental infor- 
mation concerning the turbulence and transmission of the earth 1 s atmos- 
phere and their effects on wave propagation at optical frequencies. 

They will also provide critical data needed to design future optical 
communications systems for ground to satellite applications. In 
addition a number of engineering measurements will be made for the 
purpose of evaluating the performance of the tracking and communications 
systems and confirming the feasibility of the design concepts employed. 
These measurements will include: 
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a- Measurement of the return signal strength (joules /pulse) in 
the ground based acquisition aid (GBAA) » 

b. Determination of the maximum range for acquisition with the 

GBAA. 

c . Measurement of the frequency of loss-of- track and time to 
reacquire. 

d. Observation and video-recording of a real time television 
transmission. 

e. Monitoring of certain TT house keeping data which will provide 
information on the operation and reliability of the various subsystems. 

This data will be analyzed and evaluated in conjunction with other 
flight measurements, such as scintillation, and data from peripheral 
experiments as described below. 

Engineering Analysis 

1. Objectives 

From an engineering viewpoint the high altitude aircraft 
test are intended to perform a number of functions. These include: 

a. To demonstrate the feasibility of the technical approaches 
taken to the problems of acquiring and cooperatively tracking. 

b. To demonstrate the feasibility of high data rate optical 
communication over vertical paths through the atmosphere. 

c. To demonstrate the reliability of optical tracking and com- 
munication systems under actual operational conditions not too dis- 
similar to those which would be encountered in an earth to space link. 

d. To provide engineering data which will be valuable in the 
design of future optical communications systems. 
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e. To identify the sources of error in the tracking and communi- 
cations systems and to estimate their relative importance in terms of 
overall system operation. 

In order to achieve these objectives a detailed evaluation of the 
operation of each subsystem should be performed. In the following 
section we have considered each subsystem and identified certain areas 
which we feel should be investigated and recommended measurements which 
should provide useful engineering data. In the following discussion it 
is assumed that the system is operating normally and within design 
specifications. Needless to say an important aspect of any engineering 
evaluation is the identification of those areas in which the system 
does not perform as expected or does not meet its design specifications. 
Since it may not always be possible to anticipate what these areas will 
be or what measurements should be made we have excluded measurements 
on malfunctioning systems from this discussion. Such measurements 
should more properly be considered trouble-shooting. We also exclude 
from consideration routine measurements which would fall into the 
category of acceptance-testing. 

2 . Ground Based Acquisition Aid (GBAA) 

In the acquisition mode of operation the ground based acqui- 
sition aid (GBAA) depends upon the return of a single laser pulse 
in order to detect the incoming aircraft. This makes the system very 
sensitive to the effects of atmospheric scintillation since acqui- 
sition could be missed if the critical pulse occurred at an instant of 
deep fade. The designers of the GBAA (ITT Corp.) have attempted to 
minimize this possibility by dispersing the retro-reflectors on the 
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aircraft over distances large compared to the correlation distances 
for atmospheric scintillation. In this way it becomes unlikely that 
all retro-reflectors will be in a region of fade at the same instant. 

This technique is in effect a form of aperture averaging for the 
returned signal. Considerable effort has been devoted to the analysis 
of the intensity of the returned pulses and the effect of dispersing 
the retro-reflector arrays [15]. 

It is important to the design of future acquisition and tracking 
systems that the effects of scintillation be experimentally determined 
and the effectiveness of dispersing the retro-reflectors be evaluated. 

This may be accomplished during the high altitude aircraft test by 
measuring the returned signal strength and comparing the observed data 
with the theoretical predictions. Quantities to be measured should 
include : 

a. The mean intensity of the returned laser pulse. 

b. The probability density function of the returned pulse amplitude. 

In order to compare these measurements with the theoretical pre- 
dictions both quantities should be measured simultaneously with scintil- 
lations experiments performed with the up and downlink communications 
systems. 

A second engineering measurement on the GBAA will consist of 
determining the frequency of loss-of -track and the ability of the 
system to reacquire. It is expected that due to scintillation or sys- 
tem noise the GBAA will occasionally loose track of the aircraft. When 
this occurs the GBAA will automatically switch from the track-mode to 
a reqcquire mode and the presence of this condition will be indicated 
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by a light on the GBAA control console. In order to measure the fre- 
quence of loss-of-track it will only be necessary to connect a counter 
to record the number of times the reacquire mode indicator is turned 
on. It might also be of interest to arrange a timer which would 
be started when the reacquire mode is entered and stopped when the 
track mode is entered. In this way the mean time required to reacquire 
the aircraft could be determined. 

The expected frequency of loss-of-track will depend upon the 
detection probability for a pulse and thus may be computed in a 
manner similar to that used to determine the probability of acquisition. 
Comparison with the observed frequency with which track is lost will 
therefore provide an additional verification of the atmospheric model 
used in designing the GBAA. 

3. Main Tracking System 

Errors in the main optical tracking system are expected to 
arise from three sources; atmospheric scintillation, mechanical errors 
in the telescope mount and associated drive motors, and errors in the 
control electronics. Analysis of the tracking system should not only 
determine the overall tracking accuracy of the system but also should 
attempt to identify the contribution of each of these factors to the 
total system error. 

Available data will consist of the telescope mount position encoders 
read-out, the signal to the fine pointing beam storers and the output 
signal from the tracking detector. The sum of these three signals 
will be the instantaneous bearing of the target aircraft. From this it 
will be necessary to remove the actual aircraft motion and the apparent 
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pointing fluctuations due to atmospheric scintillation. Changes in the 
measured bearing due to actual aircraft motion can be assumed to be 
slow and smoothly changing so that the variance of the fluctuations 
about an instantaneous mean can be attributed entirely to system error 
or atmospheric effects. Simultaneous observation of scintillation 
will provide a measure of the integrated atmospheric turbulence along 
the optical path from which the angle of arrival fluctuations may be 
estimated.* After correcting the variance of the observed tracking 
fluctuations for actual target motion and atmospheric effects the 
remaining variation may be attributed to electrical or mechanical 
tracking errors. While this method is at best rough it is felt that 
it should yield a fair estimate of the tracking system performance. 

4. Communications System 

The principal engineering analysis which will be performed 
on the communications system will be measurements of bit error rates 
for the downlink P.C.M. system. These measurements have been described 
in detail in the scientific measurement program document so will not 
be discussed further here. 

In addition to quantative bit error rate measurements a real time 
television picture will be transmitted to allow an evaluation of 
picture quality. By alternating between bit error rate measurement and 


& 

This estimate will not be exact since scintillation and apparent 
angle of arrival fluctuation depend upon the index of refraction 
structure constant multiplied by an appropriate lever function and 
integrated along the optical path. Since this lever function is not 
the same in both cases it will be necessary to make as assumption 
concerning the structure constant profile which may effect the 
result. 
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video transmission it will be possible to determine the bit error rate 
for each sample picture. Thus the subjective evaluation of picture 
quality amy be correlated with a measured bit error rate. Considerable 
research has been conducted to correlate bit error rates and picture 
qualities for R.F. television systems. However, because of the 
differences in the statistical nature of conventional communication 
channels and atmospheric channels it will be of interest to repeat 
these studies for an optical communications system. 

The uplink communications channel will be used only to transmit 
operational commands to the airborne experiment package. Unless it 
should prove to be unreliable, which is not expected to be the case, 
it's operation will not be well suited to any type of quantitative 
analysis or evaluation. 

Conclusions 

In attempting to document an engineering measurements program one 
encounters certain difficulties not presented by a scientific measure- 
ment program. As a scientific experiment the high altitude aircraft 
project involves certain well defined aspects of optical propagation 
which one wishes to investigate. It is therefore possible to detail in 
advance exactly what data is required and what measurement should be 
made. From an engineering point of view the most interesting and 
important aspects of the experiment may well be those in which the 
experiment does not perform exactly as expected. If a system operates 
precisely as it was designed then thete is no engineering problem of any 
importance left. This consideration makes it impossible to detail in 
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advance all the engineering experiments which will probably be con- 
ducted in connection with the project. Secondly, unlike a scientific 
experiment where much qualitative data is collected, engineering experi' 
ments involve much more qualitative or subjective observation on the 
part of the personnel involved. That is to say an engineering experi- 
ment is often more of a learning experience for the engineer than a 
formal process of data taking and therefore much less subject to 
advance planning and documentation. 

In the preceedlng sections several areas have been identified in 
which it is felt an analysis of the operation of the High Altitude 
Aircraft Experiment hardware will lead to engineering data of value in 
the design of future optical communication systems. A number of 
measurements have been suggested to obtain this data. It is expected 
that as the MSFC operating personnel gain experience with the operation 
of the system many other engineering measurements will present them- 


selves . 
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APPENDIX 5 


PATH ANALYSIS FOR AVLOC PROJECT 

The following report on the atmospheric propagation experiments 
performed during the AVLOC flight was prepared at the request of Dr. 

Joe Randall and is intended as a portion of the MSEC final report on 
project AVLOC, 

Data used in this report has been taken from "Measurement Report 
on Altitude Aircraft Test" Vol. I (AVLOC Flight 9), Vol. II (AVLOC Flight 
12) and Vol. Ill (AVLOC Flight 15), final report on NASA Contract NAS8- 
29194, by Northrop Services Inc., Huntsville, Alabama. The description 
of the data reduction procedures used by Northrop Services Inc. also 
follows closely that found in the referenced report. 

Certain paragraphs in the following report have been omitted when 
the pertinent data were not available. It is intended that the omitted 
paragraphs will be supplied by MSFC personnel who are preparing other 
sections of the final AVLOC report. 

Reference numbers are also omitted since they will be specified 
during final editing in a manner consistent with the complete AVLOC re- 
port. Likewise the equation numbers will have to be modified to agree 
with the numbering of the equations in the other parts of the report. 



PROPAGATION EXPERIMENTS 


Introduction 

The principal objectives of the AVLOC test were the engineering 
evaluation of an operational optical communications system and the scien- 
tific investigation of optical propagation over vertical paths throught 
the atmosphere. In order to fulfill the second objective extensive measure- 
ments of scintillation, bit error rates, angle of arrival fluctuations and 
atmospheric attenuation were planed. Each of these quantities was to be 
investigated as a function of receiving and transmitting aperture, range 
and zenith angle, and all measurements except bit error rates were to be 
made on both up link and down link beams. Since the AVLOC tests were to 
extend over the period of a year it was also expected that seasonal varia- 
tions in propagation characteristics and the effects of meterological 
conditions could be studied. 

The following table, reproduced from the original measurements plan, 
lists the proposed propagation experiments. 
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Table . MEASUREMENTS OUTLINE 


Quantity Measured 

1. Scintillation 
(downlink) 


2. Scintillation 
(uplink) 


3. Angle of Arrival 
Fluctuations 
(downlink) 


4. Angle of arrival 
Fluctuations 
(uplink) 


5. Bit Error Rate 


Parameters 


Analysis to Yield 


Receiving Aperture 
Range 

Zenith Angle 


Receiving Aperture 
Range 

Zenith Angle 
Transmitter Aperture 


Range 

Zenith Angle 
Receiving Aperture 


Range 

Zenith Angle 


Range 

Zenith Angle 
Beam Divergence 
Transmitter Power 


Log amplitude variance. 

Probability density function. 

Verify theoretical pre- 
dictions concerning zenith 
angle, and range dependence. 

Aperture averaging effects: 

a. Reduction of signal 
variance. 

b. Change of probability 
density function. 

Difference in uplink and 

downlink. 

Log amplitude variance. 

Probability density function. 

Verify theoretical prediction 
concerning zenith angle and 
range dependence. 

Verify no uplink aperture 
averaging . 

Effect of transmitting 
aperture size. 

Differences in uplink and 
downlink. 

Variance. 

Probability density 

Aperture averaging. 

Dependence on range 
zenith angle. 

Differences between 
and downlink. 

Variance. 

Probability density 

Dependence on range 
zenith angle. 

Differences between 
and downlink. 

Verify theoretical prediction 
of BER dependence on system 
noise and irradiance fluc- 
tuation variance. 


function. 

and 

uplink 

function. 

and 

uplink 
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Atmospheric 

Transmittance 


None 


Atmospheric transmittance 
at optical wavelengths. 


Engineering 

Measurements 


Determine ability to 
acquire and track. 
Evaluate system performance 



Early in the project's design phase it was realized that certain of 
these measurements presented greater experimental difficulties than had 
been anticipated. For example , accurate measurements of atmospheric trans- 
mission requires a precisely calibrated detector on the aircraft to monitor 
transmitted power and an exact knowledge of the beam profile. It soon be- 
came apparent that it would be very costly to implement a meaningful 
transmission measurement. Since funds were limited this experiment was 
eliminated from the measurements plan. Likewise it was decided that 
meaningful measurement of the angle of arrival fluctuations would impose 
excessive requirements on the accuracy of the tracking system and these 
experiments were also eliminated. Thus at the beginning of the operational 
phase of the project the planned measurement consisted of uplink scintilla- 
tion, downlink scintillation, bit error rates and a qualitative investigation 
of the effect of atmospheric attenuation. 

Equipment problems during the early flights severely limited the amount 
of time available for data collection and prevented the completion of 
even the reduced measurements program. Useful data was collected only on 
flights 9, 12, and 15 and during these flights the system was fully opera- 
tional only part of the time. We were therefore unable to conduct extensive 
investigation into the effect of transmitter and receiver aperture size 
nor could bit error rate measurements be made. It was also discovered 
when the data was reduced that saturation in an amplifier had destroyed 
most of downlink scintillation data. 


5 



In summary, reliable data has been obtained only for the uplink 
scintillation experiments although some qualitative information concerning 
atmospheric transmission is also available* 

Data Collected 

The following is a summary of the experimental data collected aboard 
the aircraft and at the ground station. 

Aircraft Data . The data collected aboard the aircraft is shown in 
Table 2. All data was recorded on magnetic tape. 


Table 2. Aircraft Data 


Tape Channel 

1 

2 

3 

4 

5 

6 

7 

8 
9 

10 


Quantity 

Image Disector Error, Y component 
Image Disector Error, Z component 
Beam Steerer Angle, Y component 
Beam Steerer Angle, Z component 
Vibration, Pitch 
Vibration, Yaw 
Vibration, Longitudinal 
Log Amplitude Scintillation 
Linear Amplitude Scintillation 
PAM DATA 


The image disector errors (ch. 1 and 2) are error voltages generated 
in the tracking system that are proportional to the angle of arrival of the 
Incoming beam. The beam steerer angles (ch. 3 and 4) are the voltages 
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applied to the begin steerers. Vibration data (ch. 5-7) was obtained 
from three vibration sensors mounted along mutually orthogonal axes on 
the airborne communications package* Linear amplitude scintillation (ch. 

9) is the instantaneous received power at the aircraft while log amplitude 
scintillation (ch. 8) is the same signal after passing through a logrithmic 
amplifier. Channel 10 (PAM data) consists of a number of measurements 
that were made periodically and multiplexed onto a single channel for re- 
cording. The PAM data is listed in Table 3. A complete description of 
the way in which this data was obtained can be found in the system des- 
cription elsewhere in this report. 


TABLE 3. PAM DATA 


VCO Channel 

Data 

Quantity Recorded 

3 

Beacon Presence 

ON/OFF 

4 

TV filter 

OUT /IN 

5 

Receiver Aperture 

OUT /IN 

6 

Beam Divergence 

Arc Seconds 

7 

Laser Shutter 

OUT /IN 

8 

Tracker Attenuation 

Position 

9 

Receiver Attenuation 

OUT /IN 

10 

Fine Track 

ON /OFF 

11 

60 Hz Power 

Voltage 

12 

400 Hz Power 

Voltage 

13 

28V DC unregulated 

Voltage 

14 

28V DC regulated 

Voltage 

15 

5V DC 

Voltage 

16 

15V DC 

Voltage 

17 

AOCP Temperature 

Deg. C. 

18 

Cannister Temperature 

Deg. C. 
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Ground Station Data 


This section Intentionally left blank 


DATA REDUCTION 


Aircraft Data 

Selection of Data Slices . Data was recorded on analog magnetic tape 
throughout the entire time that the RB-57 aircraft was over station. 
Oscillograms were first made from the analog tapes and Inspected to deter- 
mine those periods during which uninterrupted power was being received at 
the aircraft and the airborne tracking system was operating normally. 

Based on this inspection time slices were selected for further processing. 

Analog/Digital Conversion . Selected data was digitized and recorded 
on digital magnetic tape for computer processing. A sampling rate of 2000 
samples per second was used to provide an upper frequency cut off of 1 KHz. 
The quantization level of 1024 (10 bits /sample) was used to accomodate the 
required 60 db dynamic range. 

Probability Density and Moments . The probability density function 
(PDF) and cumulative distribution function is computed for each channel 
of the selected data slices. Based on the PDF, the first four central 
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moments, i.e., mean, variance, third moment and fourth moment are computed. 
Next, the skewness and kurtosis are computed from these moments. For normal 
(gaUssian) distribution, the kurtosis will be equal to 3 and skewness equal 
to zero. 

Test for Normality of Probability Distribution . The fact of kurtosis 
being equal to 3 and skewness being zero is not enough to establish the 
normality of any distribution. At best, it is only a quick check. The 
standard statistical test for goodness-of-fit for normality is the chi- 
Bquare test. This is a one-way test giving the upper confidence limit for 
any prescribed level of significance. One problem in using the chi-square 
test is the choice of the class number and class limits for the PDF, since 
it is known that different choices may yield different results. A procedure 
for minimizing the subjectiveness has been proposed by Mann and Wald (refs. 

, and ). This procedure computes the number of classes based on the 
number of data points used and the desired level of significance. Further, 
the corresponding class limits are chosen in such a way that the theoretical 
distribution under test in each class interval is equal. The Mann-Wald 
procedure is used in the present data processing. 

Power Spectral Density . The fast Fourier Transform (FFT) is used for 
computing the power spectral density (PSD) of the data channel of interest. 
The particular version for the PSD employed in this study is based on time 
averaging over short, modified periodograms obtained by the FFT (ref. ). 
Further , use is made of the complex— conjugate relationship of the FFT for 
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real time series, so that two real data series can be transformed simul- 
taneously in order to cut the required computation time approximately in 
half (ref. ). 

Correlation Coefficient . To determine quantitatively the degree of 
dependence between any two data series, the correlation coefficient is com- 
puted. It is defined as the normalized time-lagged product between the 
given two random time series of interest. Its value should be within 
+ 1 . 

Summary of Data Recorded . Useful data was obtained only on flights 
9, 12 and 15. The data segments that were digitized and processed are 
listed in the following tables. 
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AVLOC Flight 9 Aircraft Data 


Table 4. 

Data Piece Ho, 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 


A/D Time (seconds) 

15 

10 

50 

20 

10 

30 

40 

10 

10 

20 

50 

30 

30 

20 

50 


Total Time 


395 seconds 



Table 5. AVLOC Flight 12 Aircraft Data 
Data Piece Ho. A/D Time (Seconds) 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 


200 

70 

100 

190 

80 

80 

100 

200 

70 

190, 

80 

100 


Total 


1,460 Seconds 
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Table 6. AVLOC Flight 15 Aircraft Data 


Data Piece No. 
1 
2 

3 

4 

5 


A/D Time (Seconds) 
26 
30 
60 
30 
30 


Total Time 


176 Seconds 
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GROUND STATION DATA 


This section intentionally left blank 



DATA ANALYSIS 


Uplink Scintillation 

The uplink scintillation was measured by a tracking photodetector In 
the AOCP which measured the Instantaneous received power. The photodetector 
output was recorded in three modes; linear, logarithmic and base band. 

Ihe base band channel consisted of direct recording of the detector output. 

In order to avoid the effects of background radiation the uplink beacon 
was modulated at 10.7 KHz, and the amplitude of this frequency component 
recorded (linear channel). The 10.7 KHz signal was also passed through a 
log- amplifier before recording (logarithmic channel). Inspection of the 
data has shown that the amplitude scintillation was so large that the linear 
range of both the linear and base band channels was exceeded. The log— 
arlthemlc channel was well within the dynamic range of the system. For 
this reason we have based all analysis on the logarithmic channel data. 

The available data consist of 15 segments from flight 9, 14 segments 
from flight 12, and 3 segments from flight 15 as Indicated in tables 4 and 
5. At least one 10 second segment was selected from each data piece for 
analysis. 

A summary of the reduced data is included in appendix . 

Log Normalicy 

According to the most generally accepted theory of atmospheric tur- 
bulance amplitude scintillation is expected to follow a log normal pro- 
bability distribution, i.e. , the logarithm of the amplitude fluctuations 
should be a gaussian random variable. To test this hypotheses three statis- 
tical techniques were used, the method of moments, the chlrsquare test and plots 
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of the cumulative probability of the log and linear scintillation on 
"probability paper". In the last method the horizontal coordinate scale 
is such that if the data is Indeed normal the resulting plot will be a 
straight line. Plots for the log and linear recordings of a typical data 
segment are shown in the following figures. Other data segments yielded 
essentially similar results. 

Inspection of the probability distribution curves and the statistical 
moments Indicates that in all cases the observed scintillation more near- 
ly approximated a log-normal distribution than a gaussian. The goodness- 
of-fit to the log-normal distribution was moderately good. There was 
a general trend for the fit to become worse as the amplitude variance 
increased. We have carefully considered all possible sources of experimental 
error that might bias the observed probability distributions; a discussion 
of the possible sources of error is given in a following section. Although 
we have been unable to identify any single error that could be expected 
to cause the observed discrepancy with log— normalicy it is still possible 
that the effect might be due to the accumulation of errors from several 
sources. 

In summary, the observed data shows a strong tendency toward log- 
normalicy. On the basis of the observed data we are unable to determine 
with certainty whether the observed departure from the expected distri- 
bution is due to experimental error or represents a real deviation of the 
propagation process from the theoretically predicted 'result. 
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CUMULATIVE PROBABILITY 
DENSITY PLOT (LOG) 

Flight ^ 

Channel 8 

Data Piece 6 

\ 

Figure 1 1 

i 
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CUMMULATIVE PROBABILITY 
DENSITY PLOT (LINEAR) 

Flight 12 
Channel 9 
Data Piece 6 


Figure 2 
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Amplitude Variance 

As previously stated the base band and linear scintillation data showed 
saturation effects because the dynamic range of the signal exceeded the 
linear range of the recording channel. For this reason we have taken the 
variance from the log scintillation channel only. 

It is convenient to convert the variance of the logarithmic signal 
to intensity variance. This may be accomplished as follows. We assume 
that the log channel output is a voltage V which is related to the re- 
ceived intensity I by 

V = ctlnl + y 1 

Likewise the output of the linear channel S is given by 

s » 01 2 

Where a, 0 and y are constants that depend upon the amplifier gains in 
the two channels. Solving equation 1 for I we have 

I - exp[ V ~ ^ 3 3 

hence 

S » 0 exp[ .Y - r Y ] 4 

or 

S ■ exp + ln0] 5 

ct 


19 



we may replace the quantity -Y'/ot + InB by a new constant 6 


_ 6 . V/o 

S ■ e * e 


Now let S be the mean value of S. 
o 


S ** e^ <e^ a > 
o 


The variance of S Is 


C - <(S-S) > 

8 O 

_ .6 v/o S. v/a. v2 

Cg «= <(e e “ e <e >) > 

C s - e 2 * <e v/a > 2 D 2 > 

Now we assume that 1 la log normal » hence V Is a gausslan random variable. 
For any gausslan random variable X 

2 2 

<e ax > = exp [ a<x> + y- <(x-<x>) >J 


<e aX > ■* exp[a<x> + Cx] 


8 

9 

10 


11 

12 


a • a constant. Hence 


<e V/2 > 


exp[^~ + 

a 2 a 


26 v/a 2 
e <e > 


/ 2v/a v/o \ 


c <* s 


c - s 


\ <„ v/3 > 2 < c v/o > / 

iv v/ < . j 

l <.*'** J 
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14 


15 


16 
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then using the above relationship 



or Introducing the standard deviations^ and (T* is place of the variances 

s r 

a - and 0 " 19 

s s v v 


we finally obtain. 


But since 


oj S 

s o 



20 


and 


a 

B 
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22 

23 


which Is the required expression for the normalized standard deviation of 
the intensity fluctuations. 
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The one remaining constant in equation 20 and 23 can be evaluated from 
the calibration curve of the scintillation detector. Based on post flight 
calibration data, a has been found to be 0.532 volts. 

The following table gives the log variance and the normalized 
intensity variance Oj/I o f° r each data segment. The mean and standard 
deviation of f° r eac h flight is given in Table 8. 

The quantities shown in Table 8. are the mean and standard deviation of 
the received power and hence are proportional to the mean and standard 
deviation of irradiance. The S. D. of power is in fact equal to the S. D. 
of irradiance because of the normalization. Most theoretical results , 
however, are expressed in terms of the log amplitude variance C^. This 
quantity is easily obtained from the data as follows. is defined as 



where 1 is the log amplitude 

1 o In (A /A) = l/21n(I/l) 25 


Making use of equation 1 and the definition of Cy 


we have 



26 


27 


which upon expanding the square and collecting terms becomes 
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Table 7 

Mean and Variance for Log Scintillation (Uplink) 


Flight 

Data Piece 
No.* 

Mean 

V 

. 0 

Standard Deviation 
ov 

Normalizes Intensity 
Standard Deviation 

°i n o 

9 

2 

3.114 

.332 

.690 

9 

3 

3.139 

.325 

.673 

9 

4A 

3.185 

.320 

.660 

9 

4B 

3.206 

.329 

.683 

9 

5 

3.097 

.282 

.566 

9 

6 

3.147 

.344 

.720 

9 

7 

3.129 

.298 

.607 

9 

8A 

3.028 

.198 

.385 

9 

8B 

3.038 

.198 

.385 

9 

9 

3.307 

.312 

.641 

9 

10 

3.077 

.307 

.629 

9 

11 

3.153 

.330 

.685 

9 

12A 

3.140 

.326 ' 

.675 

9 

12B 

3.094 

.276 

.556 

9 

13A 

3.090 

.143 

.274 

9 

13B 

3.085 

.151 

.290 

9 

14 

3.271 

.302 

.617 

9 

15 

3.114 

.289 

.586 

9 

16 

3.068 

.308 

.631 

12 

1 

2.776 

.315 

.648 

12 

2 

2.689 

.434 

.972 

12 

3 


.433 

.944 

12 

4 

2.650 

.380 

.816 

12 

5 

2.743 

.401 

.874 

12 

6 

2.446 

.293 

.595 

12 

7 

3.549 

.428 

.493 

12 

8 

2.955 

.422 

.936 

12 

9 

2.966 

.334 

.695 

12 

10 

2.774 

.335 

.698 

12 

11 

0.191 

1.588 

8.605** 

12 

12 

Missing 



12 

13 

0.189 

1.773 

25.813** 

12 

14 

3.0468 

.426 

.948 

15 

1 

3.116 

.424 

.942 

15 

2 

2.976 

.479 

1.118 

15 

3 

3.021 

.499 

1.188 

15 

4 

Missing 



15 

5 

Missing 



*A and B Indicates 

two 10 sec. 

segments taken from the 

same data piece. 


**These data not used. , 
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Table 8 

Averages for Individual Flights 


Flight 9 

Average 

Standard Deviation 

Mean 

3.130 

0.070 

Normalized S.D. 

0.576 

0.134 

Flight 12 



Mean 

2.859 

0.283 

Normalized S.D. 

0.767 

0.156 

Flight 15 



Mean 

3.038 

0.058 

Normalized S.D. 

1.083 

0.103 
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likewise 


C v - o 2 ^(lnI - (ini)) 2 ') 


Cj - 1/4 ((lnl - (lnp) 2 y 


hence 


Using this expression we may calculate the log amplitude variance and 
corresponding 5. D. from the variance of the signal voltage V. The re- 
sults for each flight are given in the following table. 


Table 9 

Log Amplitude Variance 


Flight 9 
c 1 » .2873 

= .074 

Flight 12 

= .3819 
C x - .133 

Flight 15 

- .4673 
C ± » .190 


Average all Flights 
C ± » .132 


S.D. (aj) - .59 


S.D. (dj) - .051 


S.D. ( CT;l ) - .032 


25 



The observed scintillation variances are of the order of magnitude that 
one would expect from atmospheric propagation theory and an in general 
agreement with other experimental data. Minott, Bufton and Fitz Maurice 
[ ] report an average of 0.18 for data collected during GSFC Balloon 
Atmospheric Propagation Experiment (BAPE I) . This is in good agreement 
with the AVLOC value of 0.13. 

Propagation theory predicts that for an infinitely small source the 
log amplitude variance is given by [ ] 

Z 

C L (0) « .56k 7/6 f dsD n 2 (S)[Z-S)(f) J 5/6 31 

J 0 


where Z is the source to receiver distance measured along the propagation 

2 

path, k is the propagation number and Cjj (s) is the index of refraction 

structure constant which expresses the strength of turbulence at each 

point along the propagation path. A model of the turbulent atmosphere 

2 

proposed by Hufnagel and Stanley [ ] predicts that (h) will be given 

by 


S 2 (h, - s V l ».^ 


32 


where h o is a scale factor usually taken to be 1000 meters, h is attitude, 

2 

and is a constant that expresses the gross turbulence strength. Both 

2 ° 

and h Q depend strongly upon local me tieoro logical conditions. It is 

° ° 2 
generally agreed that while equation 32 predicts the gross variation of 

with altitude, the actual turbulence structure is layered and shows strong 
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local fluctuations not predicted by this model. These turbulence layers 
are associated with wind sheers In the upper atmosphere and often occur 
near the tropopause. For these reasons a proper analysis of log amplitude 
variance would require a knowledge of the turbulence profile obtained 
independently of the optical measurement. Such data could be obtained from 
radiosonde borne microthermal probes but was not available during the 
AVLOC test. 

Equation 31 must be corrected for the finite source diameter of the 
laser beam and for aperture averaging effects at the receiver. Fried 
[ ] has shown that this correction depends on the normalized source size 

ka 2 

n ■ • • — sec 0 33 

h 

o 

where is the real part of the complex gaussian beam parameter and 
0 is the zenith angle. Once 0 has been computed the required correction 
can be obtained from curves in the quoted reference. 

Fried [ ] has investigated the effects of aperture averaging and found 

that the signal variance is given by 


o 


2 

S 



pdp K o ( Pl D) C x <p) 


34 


Where D is the receiver aperture diameter, the intensity covariance 

at the receiver and K the diffraction limited optical transfer function 

o 

of the receiver. From equation 33 one may compute an averaging function 6 
which is the ratio of the observed signal variance to the signal variance 
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in the absence of aperture averaging. For an infinite plane wave pro- 
pagating from apace to earth one finds that 9 depends on the turbulence 
strength C^(0) and a normalized aperture 

D//4h sec9/, 35 

O K 

For our test condition Fried results indicate that 0 will be greater than 
0.9. This theory is not strictly applicable to the AVLOC test for these 
reasons, viz, we have a gaussian beam wave not a plane wave, the pro- 
pagation direction is from earth to space, and it neglects the effect of 
the central obscuration in the cassigrainian receiving telescope. Never 
the less Fried's results should provide an estimate of the magnitude of 
the aperture effects* Clearly this effect is small, probably less than the 
probable error in the data, and can be neglected. 

Because of the lack of precise knowledge of the turbulence profile 
that existed during the AVLOC test and the absence of a theory of aper- 
ture averaging that is applicable without extensive modification we have 
not undertaken a precise comparison of the observed log amplitude variance 

with those predicted by theory. However, using Hufnagel’s model with the 

2 

usual choice of h and , equation 31 predicts values of C^O) on the 
° o 

order of 0.1 which is in relatively good agreement with the observed log- 
amplitude variances. 

Power Spectral Density 

Power spectral density (PSD) is defined as the Fourier transform of 
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the irradiance covariance function and is equivalent to the average re- 
ceived power per unit bandwidth at each frequency . Power spectral den- 
sities have been computed for each data segment. A typical PSD is shown 
in the next figure. The PSDs from other data segments were essentially 
similar. 

All PSD' s showed a more or less linear decrease with increasing 
frequency from a few HZ to 1000 HZ , the later value being the maximum 
frequency component obtainable with the selected sampling rate. In all 
cases the power density decreased about one order of magnitude at 100 HZ 
and two orders of magnitude at 1000 HZ. 


Measurement Accuracy 

Although the observed scintillation data tended strongly toward a 
log normal distribution the lack of strict log normalicy has caused some 
concern that the data may be biased by extraneous noise. We have con- 
sidered all apparent sources of extraneous noise that might contribute 
to the observed scintillation and have concluded that none could reasonably 
be expected to account for this effect. Each of the possible sources of 
experimental error is discussed in the following paragraphs, 
a. Detector Noise 

The power signal to noise ratio S/N In the detector can be estimated 
from the photomultiplier tube parameters. S/N is given by 


S/N - 


<V 


+ ( i ^) 2 
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POWER SPECTRAL DENSITY CURVE 


Data 6 

CH 8 

FLT 12 
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y & 

Where (i ) is the mean square detector current due to the signal , (i„ ) 

B w i 

2 

that due to shot noise in the photomultiplier tube and <1^ ) that due to 
Johnson noise in the load resistor. As usual (i„ ) 2 may be neglected since 

n 2 

both the signal and shot noise terms are multiplied by the square of the 
tube gain and the Johnson noise term is not. Now 


i_ " v.. ~ ^ (2m) coswmt 
5 nv 

s 

where 

G » PM Tube Gain 

n « Quantum Efficiency 

e = electronic charge 

P a received signal power 

h » Plank's constant 

V ■ signal frequency 
0 

m d modulation lndes 
■ modulation frequency, 
and 

(ijj ) 2 - 2G 2 (i c + i d ) AV 

1 a current due to incident radiation 
c 

i d - dark current 

AV « post detection bandwidth 

i - Pet)/hv„ 
c s 
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hence 


S/N 


„ 22 
2i m 
c 


e(i +i c ) v 


where a factor of 2 has been Introduced from the time average of (i B )‘ 
we assume the following parameters: 


n - 15% 

-9 

1^ ■ 1.2 x 10 amps 
-9 

P - 6 x 10 watts 

v - 6.16 x 10 14 Hz 
8 

Av =» 5 x 10^ Hz 
m - .3 


The values of n and I d are typical of the RCA 86 44 photomultiplier tube 
used for scintillation detection in the AOCP, P is the average power mea- 
sured on all flights and m represents a worse case (minimum) modulation 
index. With these values the estimated power signal to noise ratio is 
10^ corresponding to a voltage signal to noise ratio of approximately 300. 

The shot noise induced by back ground radiation has been neglected 
in the preceeding calculations since all test were conducted at night. 

The normalized standard deviation of intensity (tfj/I ) was typically 
on the order of .5. Hence it follows that the observed fluctuation in 
the output signal due to scintillation would be approximately 150 times 
greater than those due to detector noise. 
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b. Amplifier noise, recorder characteristics, electrical pickup. 

All data channel amplifiers and recorders were designed to accommodate 
the expected data without introducing either excessive noise or biasing 
the data by reduced bandwidth or dynamic range. Throughout the flight 
test the data channel amplifier and data acquisition systems appeared to 
operate normally; we do not believe, therefore, that they represent a 
significant source of experimental error. 

The scintillation data, along with all other data channels, did shew 
a weak but persistant 56 Hz component which was been identified as elec- 
trical pickup from the aircrafts (nominally) 60 Hz power system. In our 
estimation, however, this 60 Hz pickup was not strong enough to bias 
the experimental data. 


c. Pointing 

One possible source of error is false scintillation caused by pointing 
error. Since the beam has a gausslan rather than uniform intensity pro- 
file any pointing error will cause a decrease in the irradlance at the 
receiver. If the pointing error is time dependent it will cause irradlance 
fluctuations that will be indistinguishable from scintillation. The gausslan 
beam profile may be expressed as 

t t “ e / 0 

I - I e o 
o 

where 9 is the angular beam width. Then the intensity fluctuation AI 
o 

caused by a pointing error A6 is 
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or 


. _ -A9/9 

AI » I - I - I - I e o 
o o o 


AT/T , -A9/9 

AI/I » .1 - e o 
o 


A reasonable criteria la that (AI/I ) be much less than the Intensity 

o 

variance due to scintillation. The observed intensity variances are on 
the order of 1/4 hence we might require. 

, -A9/9 nt 

1 - e o .01 


- 40/0 on 

; o ■ .99 


or 


A6/9 * .01 


That is, the pointing error should be less than 1% of the beam width to 
insure that the error in observed variance is no more than a few percent 
(about 2 1/2%). On the AVLOC test the observed pointing errors were about 
1 arc seconds and the beam width 9 q was about 180 arc seconds. We conclude 
therefore that pointing error could not cause an error of more than about 
2% in the observed intensity variance. 

The above analysis assumes a smooth, gaussian intensity profile. 

If, the beam contained "hot spots" the false scintillation Induced by 
even a small pointing error could be much worse. The possibility of this 
having occurred cannot be eliminated with complete certainty; however, 
observation of the transmitted beam has indicated that the intensity pro- 
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file is smooth and free from detectable "hot spots". 

As a final check on the effect of point accuracy we have cross 
correlated the instantaneous power received at the aircraft with beam 
steerer voltages in the ground station transmitter. If pointing was 
in fact introducing appreciable false scintillation we would expect to 
see a strong positive correlation. No significant correlation was found 


d. Aircraft Induced Turbulence 

Scintillation due to the turbulence around the aircraft is not 
significant because of the thinness of the turbulence layer and its 
proximity to the receiver. Referring to equation we see that the air' 
craft turbulence would contribute an amount. 


AC e (0) 


.56 k 


7/6 


r 


S, 2 (AC) t( Z -scf)J 5/7 


dS 


Z-AZ 


to the observed log amplitude variance. Clearly with AZ small (a few 
meters) and the weighting function [<Z-S)(f )] 5/6 approximately zero (Z-S) 
this term will be negligable. 


e. Vibration 

To ensure that the observed scintillation did not include vibration 
Induced noise the vibration of the AOCP was monitored along three orthog- 
onlal axis by piezo electro sensors. Power spectra densities of the 
vibration were computed and the vibration data was cross correlated with 
scintillation. 
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The lateral vibration (yaw) showed a strong, rather sharp resonance 
at 10 Hz and a weaker broad resonance near 100 Hz. Above 200 Hz the 
vibration fell quickly to zero. The longitudinal vibration data had a 
single strong peak at 200 Hz. Both axis also have a number of other sharp 
spikes. Especially noticable are those at 55 Hz, 110 Hz, etc. which pro- 
bably are related to the aircrafts 60 Hz (nominal) power frequency. 

None of the features in the vibration PSD could be Identified in the 
scintillation PSD; neither was their any significant correlation between 
these data. We therefore conclude that aircraft vibration has no effect 
on the recorded data. 
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